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0 (54) Title: A COMPUTER-BASED METHOD FOR MACROMOLECULAR ENGINEERING AND DESIGN 
rH 

^ (57) Abstract: The present invention relates to a system and method for engineering and designing a macromolecule. An experimen- 
^ tally determined or de novo atomic structure that corresponds to the macromolecule is identified. The atomic structure is composed 
-"^ of building blocks. When the macromolecule is a peptide or a protein, the building blocks are amino acid residues. A target subset 
^ of the building blocks in the atomic structure to be optimized is identified. The coordinates of those building blocks that are not in 
the target subset are fixed. For each building block in the target subset, a large number of potential conformers is sampled. Each 
conformer to be sampled is substituted into the atomic structure and tested against an energy function that includes the equivalent 
energy of the conformer in a reference state and comprises at least one entropic term. Combinations of conformers that best satisfy 
^ an interaction energy function are identified. 
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A COMPUTER-BASED METHOD FOR 
MACROMOLECULAR ENGINEERING AND DESIGN 

1. FIELD OF THE INVENTION 

5 The present invention relates to methods for engineering and designing molecules 

which comprise building blocks that are individually amenable to systematic variation. 
Particular areas of application include the design and development of macromolecules, for 
example, proteins, peptides* nucleic acids and polymers with desired properties such as 
stability and specificity of interaction with counterpart molecules. Specifically, the present 

10 invention relates to computer-based methods that employ search methods in the space of 
available molecules or fragments thereof which could form building blocks of a molecular 
structure and which use a three dimensional description of the structure with atomic scale 
resolution. An aim of the invention is to provide guidance to the experimental scientist who 
is not able to systematically consider all of the possible combinations of building blocks 

15 which might need to be permuted* 

2. BACKGROUND OF THE INVENTION 

Biochemistry and synthetic chemistry are replete with molecules whose structures 
consist of hundreds or thousands of atoms but whose constitution can also be thought of as 

20 that of a sequence of identifiable units, each of which comprises only a small number of 
atoms. Molecules of this sort will herein be referred to as macromolecules, a term which 
can be taken to include proteins, peptides, cyclic peptides, nucleic acids, lipids, 
carbohydrates and synthetic polymers, etc. The individual units which go to make them up 
will be called building blocks, though other terms, both generic and specific, may be found 

25 in the art. For example, the building blocks of proteins and peptides are amino acid 
residues, the building blocks of nucleic acids are nucleotides, and synthetic polymers are 
built from monomers. The term monomer can also be used in a more general sense. A 
building block, on its own, will usually differ slightly from its bound form in the 
macromolecule: the reaction with the building blocks which become its neighbors in the 

30 macromolecule structure may truncate its own structure. In this way a different term is 
often used for the free building block from its bound form. For example, amino acids are 
the building blocks of peptides and proteins whereas in their bound form they are called 
residues. The term building block, as used herein, will be understood to mean both the free 
molecule and its bound form, unless otherwise evident from the context in which it is used. 

35 The chemistry and other properties of macromolecules are often understood in terms of the 
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types of building blocks employed and the order in which they are found, i.e., their 
sequence. 

Protein or peptide engineering is a process using recombinant DNA technology or 
chemical methods to modify the amino acid sequences of natural proteins or peptides to 

5 improve or alter their function. By changing, i.e., mutating, the natural amino acid 
sequence of a protein or peptide, it is possible to alter, inter alia, its stability, substrate 
specificity, activity, and inter- and intra-molecular interactions. Changes to an amino acid 
sequence can be made on a purely random basis, or can be derived from educated guesses 
based on the atomic-resolution detail of the protein or peptide three dimensional structure 

10 provided by techniques such as X-ray crystallography, nuclear magnetic resonance, electron 
microscopy, and electron crystallography. 

The mutagenesis of proteins and peptides, even when carried out non-randomly 
using structural information, can have unexpected or undesired results. First, a mutant 
protein or peptide may not have any altered characteristics from its native counterpart 

1 5 Second, a mutant protein or peptide may acquire completely different characteristics from 
those desired. Third, a mutant protein or peptide may not be properly folded, rendering it 
unstable, insoluble, lethal, or completely non-functional. Protein engineering can thus be a 
trial-and-error process for generating a properly folded mutant protein or peptide with a 
desired function. Furthermore, mutagenesis to probe the function of proteins, so-called 

20 "site-directed mutagenesis", is a time-consuming process, involving introduction of the 
mutation into the DNA coding region, transformation of the mutated sequence into the 
appropriate cells, expression of the protein, purification, and functional assays. 

Often, desired changes in protein or peptide function, e.g., altered binding specificity 
or avidity, require the simultaneous mutagenesis of several amino acids. With twenty 

25 possible naturally occurring amino acids at each position, the number of variants that need 
to be screened is enormous. For changes at three amino acids, there are 8,000 possible 
combinations; for changes at 10 amino acids, 10 l * different amino acid sequences are 
possible. Even though the number of variants may be narrowed by making educated 
guesses based on knowledge of the protein or peptide structure, a large number of mutants 

30 may still have to be made in order to engineer a properly folded protein or peptide with the 
desired characteristics. 

Ab initio peptide and protein design presents more difficulties than the engineering 
of mutant proteins and peptides. In ab initio design, nearly all of the amino acids must be 
chosen to create a properly folded peptide or protein with a desired function, making the 

35 number of possible variants even greater than for conventional mutagenesis. Furthermore, 
if the fold of the functional protein or peptide is not well-characterized, or if the structure 
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cannot be designed based on the known structure of an homologous protein (homology 
modeling), then structural information will not be available to help narrow down those 
combinations of amino acids that are most likely to adopt the proper protein fold. 
Therefore, ab initio design of proteins and peptides by in vitro production and testing of all 

5 amino acid sequence variants is impractical, if not impossible. 

Computer-based methods for designing and engineering proteins and peptides 
should allow for the identification of amino acid sequence variants that can be 
accommodated by the three dimensional structure of the protein being mutated, thus 
decreasing the number of in vitro engineering experiments that need to be performed. The 

10 central principle is that it is far easier to consider a large number of sequence variations and 
choose the best candidates through computer simulation than it is through direct 
experimentation and synthesis in the laboratory. 

Nevertheless, computational methods are non-trivial because of the complexity of 
the problem and the quality of the primary data that is accessible for immediate use. For 

15 example, the high-resolution three dimensional structures of most small organic molecules 
are available, but those of proteins are typically at lower resolution. Further, the methods of 
modeling the weak, non-covalent forces, e.g. 9 hydrogen bonds, van der Waals interactions, 
and hydrophobic interactions, that maintain the three-dimensional structures of 
macromolecules are at present very crude. And, in general, the number of degrees of 

20 conformational freedom that are required to accurately describe the structure of a protein is 
too large to enable practical exploration of its potential energy surface. Our ability to 
reliably model small changes in a protein structure is therefore limited by several factors: 
the accuracy to which the whole structure is known; the impracticality of applying usual 
optimization methods to systems as large and complicated as proteins; and the inaccuracy of 

25 the intermolecular potential functions which are needed to model the ways in which residue 
side chains determine the three dimensional structure of the protein by aligning with one 
another. 

For these reasons, current computer-based methods for designing and engineering 
macromolecules cannot efficiently and reliably predict the accommodation of variant 
30 structures by an identified protein fold, and thus have limited utility in assessing which 
sequence variants are likely to have a desired structure and function. 

Previous computational approaches to protein engineering have been limited to 
predictions of tertiary structure from sequence, geometric rather than energetic positioning 
of side chain atoms, and prediction of favorable sites of cross-linking. 
35 In an analogous way, the exploration of nucleic acid structures is subject to the same 

complexities of mathematical modeling, as well as the combinatorial problem arising from 
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the fact that at least 4 different nucleotides can be considered for each position on die 
sequence. 

Clearly, there is a need for computer-based methods for designing and engineering 
macromolecules which utilize a more complete and accurate mathematical description of 
5 macromolecular structure than has hitherto been attempted but which employ practical and 
reasonable approximations to enable efficient execution. In this way it will be possible to 
predict the structures of macromolecular variants and enable the selection of variants having 
a desired structure and function. 

Citation of a reference herein shall not be construed as indicating that such reference 
10 is prior art to the present invention. 

3. SUMMARY OF THE INVENTION 
The present invention relates to an improved computer-based method for optimizing 
specific building blocks in the sequence set of building blocks which make up a target 

1 5 macromolecule, for example the amino acid residues of a peptide or protein. In essence, the 
central features of the invention are, given a set of substitute building blocks and a set of 
positions in the sequence of the target macromolecule, use of a plurality of conformations or 
confonners of each building block; use of a scoring function to quantify and rank the 
possible structures; use of a reference structure; and use of filtering techniques for 

20 simplifying the analysis. 

In detail, the method comprises the steps of: (a) specifying at least one substitute 
building block for each position in a set of positions; (b) determining, for each substitute 
building block, one or more candidate confonners and: (i) substituting the coordinates of 
each candidate conformer or portion thereof for the corresponding building block or portion 

25 thereof in a structure of atomic resolution of said target macromolecule; and (ii) calculating 
an intrinsic energy of each candidate conformer; (c) rejecting candidate confonners having 
an intrinsic energy above a threshold value and whose weights, computed by a partition 
function are below a threshold value; (d) calculating a pairwise interaction energy for all 
possible pairs of candidate confonners that have not been rejected in step (c) and combining 

30 the sum of the pairwise interaction energies for all pairs with the sum of the intrinsic 
energies for all candidate confonners to give a solution score; {e) determining solution 
structures, from a plurality of solution structures! which have, respectively, solution scores 
that comprise an entropic term and that are lower than a predetermined threshold solution 
score; wherein: (i) each building block in said building block set is represented in each 

35 solution in said plurality of solutions by one or more candidate confonners that each 
correspond to a substitute building block that was independently specified in accordance 

4 
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with step (a) and was not rejected in step (c); and (ii) each solution score representing a 
difference iij the summed potential energy of each candidate conformer in the solution 
structure when a candidate confoimer is substituted in an atomic-scale resolution structure 
of the target macromolecule, and when said candidate conformer is substituted into an 

5 atomic resolution macromolecular structure corresponding to a reference structure. The 
application of the method may be carried out more than once, sequentially, to obtain better 
and better solutions. The solutions may then be used as suggestions for synthetic candidates 
and those molecules which are made may then be assayed against a target of interest. 
The present invention further relates to a computer program product for use in 

10 conjunction with a computer, the computer program mechanism comprising a computer 
readable storage medium and a computer program mechanism embedded therein, the 
computer program mechanism comprising an application program configured to choose a 
set of substitute building blocks in a target macromolecule, the application program, (a) 
upon receiving a request to choose a set of building blocks for a set of positions in the 

15 sequence of the target macromolecule and, being input at least one substitute building block 
for each position in said set of positions; (b) determining, for each substitute building block, 
one or more candidate conformers and: (i) substituting the coordinates of each candidate 
conformer or portion thereof for the coordinates of the corresponding building block or 
portion thereof in a structure of atomic resolution of said target macromolecule; and (ii) 

20 calculating an intrinsic energy, of the candidate conformer, (c) rejecting candidate 
conformers having an intrinsic energy above a threshold value and whose weights, 
computed by a partition function are below a threshold value; (d) calculating a pairwise 
interaction energy, for all possible candidate conformers that have not been rejected in step 
(c) and combining the sum of the pairwise interaction energies for all pairs with the sum of 

25 the intrinsic energies for all candidate conformers to give a solution score; (e) deter minin g 
solution structures, from a plurality of solution structures, which have, respectively, solution 
scores that comprise an entropic term and that are lower than a predetermined threshold 
• solution score; wherein: (i) each building block in said building block set is represented in 
each solution in the plurality of solutions by one or more candidate conformers that each 

30 correspond to a substitute building block that was independently specified in accordance 
with step (a) and was not rejected in step (c); and (ii) each solution score representing a 
difference in the summed potential energy of each candidate conformer in said solution 
structure when a candidate conformer is substituted in an atomic-scale resolution structure 
of the target macromolecule, and when the candidate conformer is substituted into an 

35 atomic res lution macromolecular structure corresponding to a reference structure. 
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The present invention further comprises a system for choosing a set of building 
blocks in a target macromolecule comprising: a central processing unit; an input device for 
inputting requests; an output device; a memory; at least one bus connecting the central 
processing unit, and the input device; the memory storing an application program 

5 configured to choose a set of substitute building blocks in the target macromolecule, the 
application program, (a) upon receiving a request to choose a set of building blocks for a set 
of positions in the sequence of the target macromolecule and being input at least one 
substitute building block for each position in the set of positions; (b) determining, for each 
substitute building block, one or more candidate conformers and: (t) substituting the 

10 coordinates of the candidate conformer or portion thereof for the coordinates of the 

corresponding building block or portion thereof in a structure of atomic resolution of said 
target macromolecule; and (ii) calculating an intrinsic energy, of the candidate conformer; 
(c) rejecting candidate conformers having an intrinsic energy above a threshold value and 
whose weights, computed by a partition function are below a threshold value; <d) 

15 calculating a paixwise interaction energy, over all possible candidate confonners that have 
not been rejected in step (c) and combining the sum of the pairwise interaction energies for 
all pairs with the sum of the intrinsic energies for all candidate confonners to give a 
solution score; (e) determining solution structures, from a plurality of solution structures, 
which have, respectively, solution scores that comprise an entropic term and that are lower 

20 than a predetermined threshold solution score; wherein: (i) each building block in said 
building block set is represented in each solution in the plurality of solutions by one or more 
candidate confonners that each correspond to a substitute building block that was 
independently specified in accordance with step (a) and was not rejected in step (c); and (ii) 
each solution score representing a difference in the summed potential energy of each 

25 candidate conformer in said solution structure when a candidate conformer is substituted in 
an atomic-scale resolution structure of the target macromolecule, and when the candidate 
conformer is substituted into an atomic resolution macromolecular structure corresponding 
to a reference structure. 

30 4, BRIEF DESCRIPTION OF THE DRAWINGS 

Figure I : Block diagram of a computer system in accordance with the present invention. 

Figure 2: Flow chart describing the steps of data input and structure validation. 

35 Figur 3: Flow chart describing steps of intrinsic energy computation and rejection of 
retainers. 

6 
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Figure 4: Flow chart describing steps of pairwise energy computation and selection of 
solution structures in order to compute solution score. 

Figure 5: Flow chart describing steps of pairwise energy computation and selection of 
5 solution strictures by use of a genetic algorithm in order to compute solution score. 

Figure 6: Van der Waals interaction energy for two methyl (-CH,) groups as a function of 
the distance that separates them. 

10 Figure 7: Electrostatic interaction between two unit charges of equal sign, either according 
to Equations 26 and 27 and using a relative dielectric constant E, of 4.0 r„ (solid line), or 
according to Equation 28 using the same dielectric constant and a screening distance r, of 
2.0 (dashed line). 

15 Figure 8: Geometrical conditions to be satisfied in hydrogen bonding. 

Figure 9: Schematic representation of the accessible surface area (ASA) of a protein. This 
surface is defined as the surface marked by the center of a water molecule (a probe P with 
radius 1.4 A) rolling around the protein while maintaining permanent contact with the van 
20 derWaalssurfaceoftheproteinatoms. The arrows indicate that the solvation potential of 
carbon and oxygen atoms correspond to positive and negative energies, respectively, so that 
carbon atoms tend to cluster inside the protein while oxygen atoms prefer to protrude 
outside. The bold surface close to the C and O atoms arc their atomic accessible surfece 



areas. 



25 



Figure 10: Conformers generally observed around a rotatable single bond (from left to 
right: gauche trans and gauche +). X and Y represent two heavy atoms, eg., the alpha 
and delta carbons of a leucine side chain. 

30 Figure 11: Distribution of Xl (Val) and %» X, O*0 values GaUS8ian 

equations (black). Dark grey curves represent the distribution of the trans, gauche - and 
gauche + conformers. In light grey are the distributions corresponding to non-rotameric 
configurations. 

35 Figure 12: Illustration of the rotamer library concept, showing the side chain conformers of 
valine and leucine. Numbers on top of each structure are the x t (Val) and Xii Xi 
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dihedral angle values; those labelled with an asterisk were obtained during the evaluation of 
Perla. 

Figure 13: Accompanying example 1, distribution of solution scores for 1600 output 
5 sequences from Perla, applied to the SH3 domain of alpha-spectrin. Wild Type sequence is 
indicated along with score of best solution and worst solutions found. 

Figure 14. Far-UV CD spectra of the SH3 protein domains, at pH 3.5: (A) wild-type 
protein, (B) first core mutant and (C) second core mutant. (Empty circles) 278K. (Filled 
10 circles) 298K. (Filled squares) 363KL 

Figure IS. Far-UV CD spectra of the SH3 protein domains, at pH 3.5: (A to C) RT-loop, 
diverging turn and distal loop protein mutants, respectively. (Empty circles) 278K. (Filled 
circles) 298K. 

Figure 16. Urea-induced denaturation of the SH3 domain proteins, monitored through 
emission fluorescence (total fluorescence above 340nm, excitation wavelength was 280nm). 
Experimental conditions were pH3.5 and 298K. (Circles) Wild-type protein, (squares) first 
core mutant and (triangles) second core mutant. Filled and empty symbols represent the 
20 data obtained from two separate experiments. The lines are the results of a least-square 
regression fitting of the data to Equation 54. 

Figure 17. Urea-induced denaturation of the SH3 domain proteins, monitored through 
emission fluorescence (total fluorescence above 340nm, excitation wavelength was 280nm). 
25 Experimental conditions were pH3.5 and 298K. (Circles) Wild-type protein, (squares) RT- 
loop mutant and (triangles) distal loop mutant. Filled and empty symbols represent the data 
obtained from two separate experiments. The lines are the results of a least-square 
regression fitting of the data to Equation 54. 
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5. MT.T AILED DESCRIPTION HF THE INVENTION 



Section 5.1 gives an overview of the invention including a computer system for 
optiinizing a set of building blocks in a macromolecule. In Section S.2, the method of the 
present invention as implemented in Perla, the preferred embodiment of the invention, is 
35 described in brief Subsequent sections describe in more detail each step of the method of 
the present invention, with emphasis on these steps as implemented in Perla. In Secti n 
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5 3 a detailed mathematical description is given of the empirical scoring fimcti nusedto 
calculate the energy difference between an optimized conformer of a mutated target protein 
and some reference state. Section 5.4 describes a modified form of the scoring function that 
is split into template, intrinsic and pairwise terms. Section 5.5 provides a detailed 
theoretical description of the molecular mechanics potential, and of van der Waals, 

5 electrostatic, and hydrogen bonding energies, which contribute to it. Section 5.6 provides a 
detailedmathematical description of the empirical potential, calculated from changes in 
solvation and entropy of the protein chain, and which introduces an approximate desenpuon 
of the interaction of solvent with the side chain conformed. Section 5.7 describes how the 
denatured state of proteins are considered in Perla. Section 5.8 describes the generation of 

10 the rotamer library used by Perla. Section 5.9 describes energy optimization and 

elimination of incompatible amino acid conformed. Optimization routines, e.g., dead-end 
elimination and mean field theory, are detailed in Section 5.10. Section 5.1 1 describes 
re-evaluation of solvation energies, and consequently, of the scoring function, for sequences 
that remain after elimination, and Section 5.12 details output from the preferred 
embodiment of the present invention. Section 5.13 details a generalization of the method to 

15 macromolecules other than proteins and peptides. 



5.1. OVERVIEW 



20 The present invention relates to a novel method for designing and engineering 

macromolecules that utilizes an accurate and complete mathematical representation of 
macromolecular structure, in order to reliably predict how precise variants of its sequence 
can be accommodated into a desired three-dimensional (3D) structure. Herein, the 3D 
structure in question may be a specific conformation of the macromolecule itself or may be 
a complex in which the macromolecule interacts with a ligand or another macromolecule. 

25 A preferred embodiment of the present invention is known as Perla. Perla is a computer- 
based method for protein engineering that manipulates peptides and proteins in order to 
identify and sort amino acid sequences capable of folding into a desired 3D structure. 

In order to model the effect of a small change on a protein structure, it is not always 
necessary to reanalyze its entire structure. Consequently, well chosen detailed structural 
information about the site of mutation may be employed to focus attention on the area of 
interest. Structural flexibility of a protein may be thought of as a large-scale consequence 
of the conformational flexibility of the building blocks of which it is composed. Here we 
may exploit the fact that residue mutation in a protein is effectively a side chain substitution 
which leaves the backbone unperturbed. That is, the part of the structure of the ammo acid 
35 building block which changes from one t another is the side chain alone. Correspondingly. 
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coiiformational analysis may be simplified by using sets of known favorable sidecbain 
conformations instead of carrying out an unconstrained energy minimization. 

To identify sequences that are most able to fold into the target protein 3D structure, 
an energy or scoring function is established, referring to a conformational state representing 
5 either the denatured or unfolded protein. A random coil or extended structure can also be 
used as reference. The scoring function, as a difference between an energy of the protein 
folded structure, and an energy of the protein denatured or unfolded state is correlated to the 
stability of the folded protein structure. 

Although the invention is described herein below mainly with application to 
10 proteins and peptides, as will be evident to a skilled artisan, the teaching herein can be 
readily adapted for use with other raacromolecules such as nucleic acids, carbohydrates and 
other polymers. 



15 5. u . A SYSTEM FOR OPTIMIZING STRUCTURAL 

TTNTTS OF A MAf!ROMO L.ECULE 

The present invention addresses the ability to engineer derivatives of large 
molecules through systematic variations of their structural components by presenting 

20 solution scores of preferred solution structures after rigorous analysis of a user-defined 
search space. The invention, as shown in figure 1, comprises a system 100 for optimizing a 
set of structural units in a target macromolecule, comprising at least one central processing 
unit 102, a user interface 104 for inputting requests, an output device 106, a section of main 
memory 1 10. and a bus 108 connecting the central processing unit, the memory, the input 
device, and the output device. Optionally, the system is connected to a network interface 

25 114, via which access to a Protein Database 1 16, such as the PDB, can be achieved. The 
memory stores an operating system, 120, a file system 122, at least one set of molecular 
mechanics parameters 124, a cache 126 and an Application program 130, also known as 
Perla. In a preferred embodiment, system 100 comprises an multi-processor computer, 
capable of carrying out calculations in parallel. 

30 Application program 130 is configured to optimize a set of structural units in the 

target macromolecule and comprises a number of modules, including but not limited to: a 
solution score module 132, to compute a solution score of a mutated or target structure; a 
sub-rotamer module 134 to calculate distributions of sub-rotamers for rotamer side chains; 
an input reader 136; optionally, a library of pre-computed rotamer conformations 138; a 

35 solvati n energy module 140 including functions that assign residues to be "core" and "non- 
core"; a surface area module 142 for computing the solvent accessible area or other surface 
of a target structure or mutated form thereof; a dihedral angle module 144 for c mputing 

10 
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backbone dihedral angles on the target structure or mutated form thereof; a penalty function 
module 146 for computing a statistically derived penalty function; a mean field theory 
module 148'for optimizing side chain distributions and for computing side chain 
contributions to the entropy, a genetic algorithm module 150 for selecting an optimum 
family of mutated sequences and/or conformers; and a dead end elimination module 152 for 
eliminating unfavorable sequences and/or conformers; and a collection of protein 
fragments, 154. Additional modules, not shown in Figure I may comprise a module to 
assign atom types to atoms in an input structure, in order to assign molecular mechanics 
force field parameters and modules containing geometric manipulation utilities, such as 
may be required to translate rotamer coordinates from a library to the template. 

The macromolecule. which may be a peptide, protein, strand of DNA or RNA or a 
carbohydrate, or any organic molecule which consists of identifiable distinct structural 
units, has a 3D structure whose geometry is known to atomic resolution. The Application 
Program 130, upon receiving a request to choose a set of substitute building blocks for at 
least one set of positions, utilizes, for each substitute building block, one or more candidate 
1 5 conformers. For each detennined candidate conformer. the application program substitutes 
a building block at a position in the target macromolecule with the candidate conformer and 
calculates an intrinsic energy term for the candidate conformers. The application program 
subsequently rejects candidate conformers having an intrinsic energy above a threshold 
value and according to whether the statistical weight of the conformer, calculated from a 
7Q partition function, is below a threshold value. The application program calculates a 

pairwise interaction energy term for all possible conformers that have not been rejected by 
the threshold value criteria. The method enables determination of solution structures, that 
are ranked by solution score. In some embodiments, for example those that utilize dead end 
elimination, the solutions include a best solution corresponding to a global minimum energy 
conformation. Each building block in the set to be optimized is represented in each solution 
25 by one or more candidate conformers that were not rejected by the threshold criteria when 
substituted into the strocture of the target macromolecule. Each solution score can be 
expressed as difference between the summed potential energy of each candidate conformer 
substituted into the target structure and the same conformer substituted into a reference 
structure. The solution score comprises molecular mechanics energy terms (van der Waals, 
hydrogen bonding and electrostatic) and terms corresponding to an empirical potential 
(entropy and solvation) along with a user-defined statistical term. 

This system, when operated in a laboratory environment can provide an efficient and 
useful method of directing experimental efforts towards engineering sequence variations in 
a target macromolecule. Said system, being capable of quantifying the potency of a 
35 plurality of sequences and thereby selecting a small number which would be worthy of 
synthesis, can operate in tandem with experiment to optimize properties of interest f the 
target macrom lecule. 

11 
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5.1.2. THE PROTE IN STRUCTURE 

The computer-based method of the present invention uses an "inverse folding" 
approach, Le., aprotein backbone is chosen a priori as the native state to be designed and is 
kept fixed throughout the calculation. Fixed protein backbone atoms are the central carbon 

5 atom, Co, and the amide group, C(=0)NH, of each residue. There is no restriction that me 
protein should consist of a single contiguous backbone; Perla can accept multiple 
backbones as input. The choice of a protein topology depends on the application of the 
engineered protein. Due to the absence of backbone motions during the evaluation of 
protein sequences, it is preferable for the main chain target conformation to be correctly 

. 0 constructed from the start. In a preferred embodiment, a protein with a well characterized 
protein fold or high resolution three dimensional structure is chosen {e.g., from amongst 
those found in the Protein Data Bank (PDB), available from Research Collaboratory for 
Structural Bioinformatics(RCSB). web site address hrto://www.rcsb.org/pdb_/). As used 
herein, the "resolution" of a macromolecular three-dimensional structure is the minimum 
separation two atoms can have and still appear to be distinct and separate. Thus, the higher 

15 the resolution, Le. the smaller the separation distance at which two atoms can be 

distinguished, the more accurately determined is the structure. In a preferred embodiment, 
the protein model is solved at atom level resolution around the site of interest and the fixed 
backbone has been refined to euminate steric clashes and unfavorable main chain dihedral 
angles. For this purpose the structure may have been obtained from X-ray crystallography 
or from NMR studies. In general, the parameters employed by the user of the invention 
may be chosen to best suit the quality of the data. In a second embodiment, e.g., de novo 
protein design, the protein structure is not available or the protein fold is not well 
characterized, and methods for the construction of novel protein backbones are employed 
(e.g., WHATJF; Vriend, 1990. J. Mot. Graphics 8:52-57; INSIGHT; Abagyan el al., 1994, 
J. Camp. Chem. 15:488-506). 

25 

The 3D structures of other macromolecules can similarly be obtained from X-ray 
crystallography or may themselves be the outcome of mathematical or computational 
simulation. 



5.1.3. THE AMINO ACID SIDE CHAINS 

The computer-based method of the present invention uses a three-dimensional 
atomic description of the system to be engineered. The main chain atomic configuration 
being provided, the method is used to reconstruct amino acid side chains. The side chains 
35 of the twenty naturally occurring amino acids are b und to the backbone Cce atoms. 
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A custom-made library of discrete side chain conformations ("retainers") for each 
amino acid, compiled using dihedral angle (x„ Xv X>. &) data {rom available structures 
(preferably from those deposited in the PDB). is employed by the method of the present 
invention. The library of amino acid side chain conformations is preferably made by fitting 
occurrences of side chain dihedral angles for each amino acid side chain in known protein 
5 structures to Gaussian distributions. Since stereochemical rules were not used to generate 
the library, it contains side chain conformations that are not very abundant but are 
nonetheless important components of protein structure. Furthermore, because each dihedral 
angle is described by a Gaussian distribution, the observed range of oscillation of each 
angle is also incorporated into the library. 

1 ° In application to polymeric structures other than proteins, it may be possible to 

derive confonner libraries from means other than by direct comparison with crystal 
structures. For example, stereochemical rules may be adequate for the hydroxyl groups of 
sugar molecules; computer simulation may be most appropriate for the modeling of 
nucleotide conformations. In some circumstances, the building blocks may have 

15 msufEcient conformational flexibility to demand construction of confonner libraries. In 
such cases, the application of the method is a lot more straightforward than described 
herein, there being fewer conformers per building block. 



20 5.1.4. SELECTING POSSIBLE SE QUENCES 

The computer-based method of the present invention executes successive trials to 
consider the immense variety of sequences that can be generated as a result of protein 
mutagenesis, Le., substitution of one amino acid side chain with a different amino acid side 
chain at a given site in the protein. 

25 

In one embodiment, the user specifies which residues in the protein are to be altered. 
These may be specified by position in the sequence of the protein. To achieve this me user 
may employ specific knowledge about the 3D structure of the protein. For example, the 
user may choose residues which seem to be critical to folding or which are important in the 
definition of a binding site. In a preferred embodiment, Perla itself, through use of its own 
30 scoring function (see below) may automatically identify the building blocks which are to be 
varied. In either case, it is not necessary that the selected residues form a contiguous stretch 
of the sequence of the target macroraolecule; nor is it necessary that any pair of the selected 
residues is adjacent in the sequence. 

35 In one embodiment, the user may also specify a list of possible mutations i.e., 

residues to be considered at each residue position in the protein or a broad category of 
desirable mutations, m another embodiment, Perla may analyze the immediate 

13 
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environment of the selected building blocks and choose mutations which arc likely to cause 
the least disruption to that locale. For example, it may be appropriate to consider only 
"polar" amino acids at a particular position which is already occupied by a polar sidechain. 

Sequence sampling as embodied in the method of the present invention consists of 
searching the required amino acid side chains within the rotamer library and fitting these 
onto the chosen backbone. Side chains of the amino acid residues that are not mutated can 
remain structurally fixed or be moved, as desired by the user of the method. 



10 



5.1.5. SPORING THF SOLUTION STRUCTURES 



In order to evaluate the degree of fit of a combination of amino acid side chain 
rotamers to a protein structure, the method of the present invention utilizes a scoring 
function made up of a sum of terms which gives rise to a solution score. Unlike previous 
methods for protein modeling, the method of the present invention not only considers the 
15 global sum of these terms, but also requires that individual terms satisfy constraints found 
in natural proteins. 

Because the nature of the application of the method is to produce a number of 
different structures, each of which is distinct from the target, it is not possible to 
meaningfully compare the energies of each Consequently, the use of a reference structure 
20 for each separate solution structure enables the structures to be compared in terms of their 
inherent stabilities. In a preferred embodiment, in order to assign a solution score to a 
particular solution structure, the method calculates the difference in potential energy 
between the mutated protein structure and a reference structure whose sequence is the same. 

A preferred embodiment of the present invention calculates a potential energy for 
the native and denatured (reference) states. For the latter, in a preferred embodiment, 
sample conformations are taken from structures present in the PDB; this method is 
described in more detail later. The energy difference between the two states serves as a 
score, and the higher the score, i.e., the larger the energy difference between the two states, 
the better the degree of fit of the chosen sequence to the overall native-state protein 
30 structure. 

The estimation of the native state potential energy requires that the optimal 
association of amino acid rotamers be found. For peptides longer than a few residues, an 
exhaustive sampling of every possible combination of rotamers is not practical. Choosing 
the most likely organization of side chains is a significant combinatorial problem and, 
35 therefore, the method of the present invention employs optimization routines. The 
underlying principle of available optimization methods, e.g., dead-end dirnination and 

14 
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mean field theory, is that the energy is expressible as a scoring function comprising a term 
to describe the fixed template, one sum of terms intrinsic to every single amino acid of the 
sequence and a second sum for all pairs of residues. 

In the preferred embodiment of the present invention, a user-defined set of rotatable 

5 side chains is modeled in the context of a fixed collection of atoms, which include main 
chain atoms and the side chain atoms of residues that are not included in the modeling set. 
Together, the fixed atoms are the template, the structure which is the direct environment of 
the side chains that are subject to modeling. The calculation of the sequence-independent, 
constant energy term corresponding to the template is not required for the evaluation of the 

10 optimal set of side chains, but can be determined in order to estimate the quality of the 
template structure itself. Both the intrinsic and pairwise energy terms are similar in nature 
and are established to correlate with observed structural parameters in proteins. The 
intrinsic energy term arises from interactions between the (fixed) template and the 
(rotatable) side chains, while the pairwise energy term arises from interactions of the 
(rotatable) side chains amongst themselves. Additionally, both the intrinsic and pairwise 

15 energy terms contain contributions which depend only on the nature of the residue. A van 
der Waals component associated with the packing of atoms, an electrostatic tenn associated 
with ion pairs, and a hydrogen bonding term, are contained in both the intrinsic energy term 
and the pairwise energy term of the scoring function. In the preferred embodiment of the 
present invention, not only the global sum of the scoring function that is considered, but 

20 also each individual term must satisfy determined constraints, as reflected in naturally 
occurring protein structures. 

In other embodiments of the present invention and in application to macromolecules 
other than proteins, it may be preferable to use more than one reference state for the 
calculation of the scoring function. Alternatively, in other embodiments, the use of a 
25 reference state may be unnecessary. 

In yet other embodiments, the scoring function may comprise a term quantifying the 
interaction between the macromolecule and some binding partner. For example, the 
macromolecule may be an enzyme and the partner its substrate; in another example, the 
macromolecule may be a peptide sequence and its partner may be another peptide sequence; 
30 in a further example, the macromolecule may be a nucleic acid and its partner may be a 
protein or some fragment thereof. 



5.1.6 OPTIMIZING OR EVALUATING THE S OLUTIONS 

3 The combinatorial problem of side chain building is solved in the method f the 

present invention by calculation of mean field energies. This novel integration of Mean 

15 



BNSDOC1D: <WO 0116810A2_I_> 



WO 01/1681 



PCT/EPOO/08S04 



Field Theory, an iterative approach, into a protein modeling method provides a measure of 
the entropy f the molecule and allows for the c nsideration of all possible amino acid side 
chain confonners rather than just the global energy minimum, which is a more accurate 
description of macromolecular structure. 

5 The foregoing steps are applicable to each individual substituted residue; the 

problem of considering many alternative substitute residues at many different sites is also 
combinatorial in scope. The invention addresses this aspect by the technique of "dead end 
elimination*' in which certain candidate rotamers may be eliminated from the search space if 
their energy scores obey certain inequalities with respect to the scores of the other rotamers 

10 present in the same solution. Consequently the overall invention comprises two distinct 
methods of addressing problems of a combinatorial complexity. 



52 TUT STEPS OF SEQUENCE M ODELING AND EVALUATION 

USING PERLA 



15 



As described in Figures 2, 3 and 4, Perla, the preferred embodiment of the present 
invention, first reads the user-specified input at step 228. The input comprises at least four 
pieces of information. The first is the protein structure 222, i.e.. the atoms comprising the 

20 specified template (or "target") protein conformation and foeir Cartesian coordinates. These 
coordinates may have originated as fractional coordinates from a Protein Data Bank (PDB) 
file. There is no restriction that the atoms comprising the template form a connected unit, 
i.e., the protein may have multiple backbones, or several discrete proteins or peptides in 
juxtaposition may constitute the template. It is preferred to utilize other computational tools 

25 which ascertain the appropriate protonation state of the residues and ionize them as 

applicable, before passing the coordinates to Perla. A second item of user-specified input is 
a selection of amino acids 2 14 to engineer. A third item of user-specified input is a list of 
positions 210, or an indication that Perla should make some determination of its own in mis 
regard. In one embodiment, the user can specify that the side chains of just those residues 

30 in ihe set are subject to optimization. In another embodiment, Perla can automatically . 
determine which residues in the vicinity of those specified should also be subject to 
optimization of side chain conformations. In another embodiment, Perla can optimize all 
side chains. A fourth item of user input is a series of adjustable input parameters 226 that 
set weights for the different energy terms, place thresholds and penalties to control the flow 

35 of output and tune the effectiveness of the optimization procedure. 
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With this information) at step 230, Perla places amino acid residues from list 214 
into the protein structure at positions specified in the list of residues 210. Side chains that 
correspond to the list of amino acids to model are obtained from a rotamer library 232. 
This structure is an initial form of the mutated structure and it can be validated, at step 234, 
5 by calculating various terms in the solution score. 

As described in Figure 3 the intrinsic energy term of the scoring function, i.e., the 
intrinsic energy of interaction between the rotamers and the template protein structure, is 
computed, for validated structure 300 step 302. The intrinsic energy is determined by 
summing and optimizing van der Waals, electrostatic, and hydrogen bonding energies, 
computed utilizing molecular mechanics force field parameters 308 that may be read from a 
separate file or specified by the user. Those skilled in the art will recognize that, although 
the set of parameters used by the method of the present invention relate to amino acids, 
corresponding parameters for nucleic acids and other organic compounds, e.g., 
carbohydrates, are available and can readily be integrated into the method as applied to 
other types of macromolecules. 

Main chain entropy costs, vibrational side chain entropies, intrinsic contribution to 
the solvation energies and a penalty function are added to the molecular mechanics terms. 
Intrinsic energies for both the mutated structure 312 and a reference structure 306 are 
computed. In a preferred embodiment, the intrinsic energies for each side-chain and all of 
the contributions to the intrinsic energy, including, for example the solvation energy for 
each side-chain, are stored in memory. 

In step 3 14, one or more selection criteria are applied to side chain confoimers and 
those that are not compatible with the template structure are abandoned 310. In a preferred 
embodiment, a first criterion is an energy threshold, above which rotamers are rejected and 
a second criterion is a weighting from the partition coefficient, below which rotamers are 
rejected. 

Subsequently, as described in Figure 4, at step 402 pairs of selected rotamers 400 
that were not rejected at step 314 are considered in order to evaluate the pairwise (side 
chain-side chain) component of the scoring function. As for the intrinsic energy 
computation, molecular mechanics force field parameters 410 are utilized in order to 
compute van der Waals, electrostatic, and hydrogen bonding energies. Molecular 
mechanics terms are summed and optimized, 6tep 404, and then a pairwise solvation 
contribution and vibrational entropy and statistical penalty is added for both mutated 
structure 412 and reference structure 406. No elimination need be performed at this step, 
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since the identification of an energetically disfavored pair does not necessarily imply that 
the participating side chains are incompatible with the target protein fold. In a preferred 
embodiment, the pairwise energies for each pair f sidechains and all of the contributions to 
the pairwise energy for each pair of side chains, including the solvation energy contribution 

5 for each side-chain, are stored in memory. Although this can be a very large matrix, it is 
better to try to store it than to repeatedly recompute the pairwise energies. For a very large 
contribution, for example comprising many variable side chains and utilizing a genetic 
algorithm, an upper limit threshold is placed on the number that are stored. Li one 
embodiment, this may be 1 Gigabyte of random access memory. In a preferred 

10 embodiment, when a multi-processor machine is available to run Perla, calculation of the 
pairwise energies can be carried out in parallel. 

In order to reduce the number of sequences to sample and to ultimately find 
that which has the optimal sequence-to-structure relationship, e.g., the lowest native state 
potential energy (lowest score) or the greatest energy difference in energy from the 
15 reference state, dead-end elimination is used to establish which sequences cannot achieve 
the energy minimum, step 408. Those sequences that cannot achieve a desired energy 
minimum are marked and abandoned, 414. In an alternative embodiment, step 408 
comprises the selection of a subset of sequences instead of a rejection and is achieved with 
the use of a genetic algorithm which gives a predefined number of best sequences. In a 
preferred embodiment, when a multi-processor machine is available, the dead end 
elimination can be run in parallel. 

In one embodiment, the intrinsic and pairwise energies can be added to one another to 
give an estimate of the stability of the structure. This is mainly of interest for small 
25 macromolecules or where only a small number, say up to 3 residues, are being mutated and it 
is not worth the overhead to go through a dead end elimination, genetic algorithm or mean field 
theory calculation. 

For all remaining sequences, mean field theory, step 416, enables the estimation of 
30 weights for all side chain retainers, which are then used to compute the solution score 41 8 of 
each sequence. Sequences that do not score well are rejected 420, while for others, the 
solvation term is re-evaluated 422. Some sequences may also be eliminated at this step if they 
have poor solvation energies. 

Finally, in the output 424 from the program, the solution structures of the engineered 
35 sequences are accompanied by a description of the energy terms that contribute to the scoring 



20 
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iluiction and a set of three-dimensional Cartesian coordinates that describe the modeled 
solution structure. 

In two further embodiments, shown in Figure 5, a genetic algorithm is utilized. In both 
embodiments, the pairwise energies are not pre-computed, but are built up 512 through.the use 
of the genetic algorithm. The selected rotamcrs 400 are passed directly into the genetic 
algorithm, 510, which simultaneously optimizes the pairwise energy 502 and the sequences. 
The usual molecular mechanics force field parameters 504 are used. Those contributions to the 
pairwise energy that have not been computed when the Mean Field Theory calculation 5 16 is 
to be carried out, are then computed 514, so that a complete set of pairwise energies is stored. 
In one embodiment, one hundred sequences from the genetic algorithm for example, are 
subjected to Mean Field Theory. The difference between the two further embodiments shown 
in Figure 5 that utilize a genetic algorithm is principally in the way in which the solvation 
energy is treated. In one embodiment in Figure 5, the intrinsic and pairwise contributions to 
the solvation energy are utilized within the genetic algorithm and subsequently in the Mean 
Field Theory calculation. Once a set of sequences has been derived, the solution scores 518 
are optionally presented. In a preferred embodiment, the intrinsic and pairwise contributions 
to the solvation energy are subtracted out and the final solution scores computed with refined 
solvation energies. In an alternative embodiment, the intrinsic and pairwise energy 
contributions to the solvation energy are not utilized in the genetic algorithm calculation. 
Instead, the refined values of the solvation energy are utilized throughout. In a preferred 
embodiment, when a multi-processor machine is available to run Perla, the genetic algorithm 
can be run in parallel. 

In yet another embodiment, a genetic algorithm may be applied to those sequences that 
have been obtained through the use of a dead end elimination step. In a further embodiment, 
it is possible to utilize a genetic algorithm to derive families of sequences without separating 
out the intrinsic and pairwise interaction energies. In such an embodiment, a genetic algorithm 
utilizes single composite solution score, to evaluate the molecular mechanics interactions, the 
solvation energies, and the statistical penalty. 

30 



15 



20 



S3. THE GENERAL FORM OF THE SCORING FUNCTION 

Central to the operation of Perla is the use of an scoring function to calculate the energy 
35 difference between a mutated version of the target protein and some corresponding reference 
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state. The way in which the reference state is constructed for proteins is described in more 
detail below, section 5.7. 



5 The contributions are regarded to be either components of a molecular mechanics model 

or part of an empirical description of solvation and cntropic factors. The theory behind each 
of these categories is described later in this section. 

In one embodiment, the general form of the solution score consists of 6 contributions, 
10 as shown in equation (1): 

P , _ ^Molecular Mechanics . E-Entrooy 

£ Scoring Function ~ & r Maintain 

Vibrational Entropy jrEntropy . m 

+£side Chain „ . , + 1 *S«»a»' n W 
+ £Solvation + ^a*}^* 1 Penalt y 

15 



Each of the components in equation (1) includes as a coefficient, a user-defined 
weighting factor, w, which can be adjusted to suit different applications. The reason for this 
20 is that some of the terms overlap with one another in the contributions that they model, and 
therefore represent over-estimates of the contributions. In a preferred embodiment, w for 
solvation is 1 .0 and all other w's are set to 0.5. Depending upon the application of the program, 
the coefficients can be adjusted to achieve a desired result. The explicit form of the terms is 
as follows. 

25 

The molecular mechanics term describes long-range interactions between pairs of atoms 
and supplies the difference in such terms between the target protein and the reference structure. 
In one embodiment, the molecular mechanics terms comprise three types of contribution, van 
der Waals, electrostatics and hydrogen bonding terms, each ofwhich is multiplied by a separate 
30 coefficient w**, w* and yf^, respectively, see equations (2a,b,c). These coefficients may be 
adjusted to permit force field parameters from different sources to be used. For example, some 
published force-fields do not have separate hydrogen bonding terms. If parameters from such 
force fields are used, w hb can be set to 0. 
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In practice the molecular mechanics terms in equations (2a,b,c) are averaged over all 
the retainers of each amino acid residue. 



The second term in equation (1 ) describes the entropy cost of fixing the main chain at 
a physical temperature, r phyi , as shown in equation (3). The temperature is typically in the 
35 range 278-328K, and in a preferred embodiment is 298K. This term is akin to a secondary 
structure propensity term and is computed as part of the intrinsic energy contribution. For each 
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residue, i, in expression (3), the contributions represent effective averages over rotamers 
because they depend only on the identity of the amino acid. 



^ amino acid/ 



-^mShain^^phys X ^ 
all residues i 



subspaccs 20°x20- 
close to I ' 



amino acid i 



(3) 



10 

The form of expression (3) is chosen so that it resembles an effective free energy term, 
AG, given by -RT InK. The third term in equation (1 ) represents the entropy cost of restricting 
the "vibrational" freedom of rotamers and is shown in equation (4). It allows priority to be 
15 givento side chain rotamers that can freelyrotate within a space corresponding to the Gaussian 
distributions determined during the creation of the totamer library. The w, are obtained from 
a partition function over the sub-rotamcrs of the rotamers. The w, are obtained from Mean 
Field Theory. The vibrational entropy term is calculated according to expression (4). 
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The fourth term in equation (1) represents the entropy cost of placing amino actd side 
chains bto the template structure where they are more hindered due to the compact protein 
environment This term, expressed in equation (5), is again a difference between the enttopiea 
of the side chains in the template and the reference. The values of the w, are obtained from 
Mean Field Theory. 



35 
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Hie fifth term in equation (1) is the solvation energy. The solvation energy is computed 
as a difference in the energy of interaction between the target structure and surrounding solvent 
and the reference structure and surrounding solvent, as shown in equation (6). 
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In a preferred embodiment, the solvation energy is computed over only one 
conformation of the solution structure, typically that obtained by using the most favorable 
retainer of each mutated residue. 

25 The last term in equation (1) is a statistical penalty function which is related to the 

identity of the amino acid residues in the sequence, as shown in equation (7), This term is 
introduced to drive the sequence design towards a sequence subspace known a priori to be 
plausible. It is not computed with respect to a reference structure of the mutated protein. 
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(7) 



We note that the entropy of the side chains is dependent upon the weight distribution 
calculated by the mean field approximation routine (section 5.10), Hence, that part of the 
35 energy is not included at all in either the intrinsic or pairwise descripti n. By contrast, the 
vibration entropy cost is used to penalize rotamers whose interaction energy (either intrinsic 
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10 



or pairwise) is only optimal for a few of the sub-rotamer conformations they can adopt (see 
sections. 8). 

For example, the values, /£fno acid / • of ^ term can represent amino acid relative 

abundance in the protein database (PDB) or can be obtained from sequence alignments related 
to the family of proteins containing the target protein. The effective temperature, 
associated with this term, might differ from the actual physical temperature 7^ used for 
entropy related terms. The factor RT^ has been introduced to equation (7) to ensure that the 
penalty term as a whole has dimensions of energy. In another embodiment, equation (7) is 
extended to include terms that describe the relative occurrence of pairs of amino acids. For 
example, during sequence alignment, it may be found that 1 0% of the sequences have alanine 
and valine at a pair of positions, but the natural abundance of alanine and valine, respectively, 
may be each less than 10%. An appropriate factor, P t for a pair of residues would be the 
abundance of the pair (in this case 0. 1 0) divided by the product of their individual abundances. 
In this way, where two residues are especially favored at a pair of positions, the value of P is 
greater than 1.0 and its log is a positive number. Accordingly the statistical penalty will 
contribute negatively, i.e., favorably to the energy. Where the value of P is 1.0, i.e., the pair 
of residues is neither favored/nor disfavored, the penalty does not contribute at all because the 
log of 1 .0 is zero. Pairs of residues found to be a fixed distance apart in a sequence are often 
found to be close in contact in the folded structure. In yet another embodiment, equation (7) 
can be further extended to include to include tenns describing the relative occurrence of triplets, 
quadruplets, etc., of amino acids. 

In another embodiment, for optimization with a genetic algorithm, an additional 
25 statistical penalty may be used, equation (8). 



15 
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30 Where N PDB is a probability of occurrence in the PDB, and N is the probability of occurrence 
in the sequence of the mutated sequence. The constants K, are user-defined weights, which are 
adjusted to give a contribution in the range of IKcal/mol' 1 . 



5.4 MODIFI ED FORM OF THE SCORING FUNCTION 
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Additionally, for the purpose of selecting, rejecting as optimizing conformers, various 
components of the scoring function are partitioned into "template", •intrinsic" and ^pairwise" 
terms. 

The context in which the scoring function should be viewed is that the protein 
comprises a fixed backbone (or backbones) of amino acid residues, some specified subset of 
which are to be varied. The backbone and the constant residues (including their side chains) 
together form the template. The side chains of the variable residues interact with the template, 
giving rise to the "intrinsic" energy term and amongst themselves, giving rise to the "pairwise" 
energy term. In the preferred embodiment of this invention, the scoring function is therefore 
decomposed into a sum of terms, described respectively as "template", "intrinsic" and 
"pairwise" equation (9). Each of these terms partitions into summed contributions equation 
(10). 
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^sequence-to-structure ^template ^intrinsic + ^painvise C) 

^sequence-to-structure ~ -^template + S -^intrinsic + ]C ^pairwise (10) 

all residues, i all residues, i> j 

In equations (9) and (10), Sequence- to- structure is an energy that is computed during 
the steps of optimization, prior to evaluation of the solution score for the mutated structure. 

In a preferred embodiment, amino acid residues are classified as "core" or '"non-core" 
prior to computing the scoring function. This classification is relevant for computation of the 
electrostatic interaction energy, through use of a variable dielectric constant, and in the 
calculation of the solvation energy. 



In other embodiments of the present invention, the scoring function may be partitioned 
into terms corresponding to interactions between atoms on the macromolecule and atoms on 
30 some binding partner of interest. Such terms may be printed in the output 



5.4.2. THE TEMPLATE TERM 
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In a preferred embodiment, in which all atoms of the template are fixed, only the 
molecular mechanics, main chain entropy and statistical penalty terms are computed, as shown 
in equation (11). 



t _ in Molecular Mechanics ^ Main Chain Entropy _ 7-^StatisticalPenalty n 

' Template ~ & Template & Template ~ & Template V 1 U 



In equation (11), the subscript 'Template" means all atoms of the template. In an 
embodiment in which all atoms of the template are fixed, the entropy contributions for the side- 
chains of the residues are not determined. The template term does not find application in 
sequence prediction where the sequence of the template is the same for all mutated structures. 

The template term is useful in structure validation and because when added to the 
intrinsic and pairwise terms, gives an Energy Contribution that is proportional to the size of the 
protein being mutated. 



20 

5.4.2. THE INTRINSIC TERM 

The invention provides an intrinsic energy term for all candidate retainers, that 
represents the interaction of each with the main chain (and any other side chain that is kept 
fixed). The intrinsic energy can be used to reject unfavorable side chain ro tamers. 

25 

The intrinsic energy term consists of S contributions, each pertaining to interactions of 
the side chains of the variable residues with the template. 
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is T ^Molecular Mechanics + E^°*Z - (12) 

C*rin- Main Oah ^MainOuun ^ 

y-iVibfiaional Entropy , *-i Solvation ^Statistical Penalty 
+ & Side Chain + ^ Side Chain* ^ Residues 

35 

The terms in equation (12) mirror those in the general expression for the energy* equation (1)- 
The intrinsic energy is expressed as a summation over all mutated residues and all those that 
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are flexible. The molecular mechanics contribution to the intrinsic energy for a candidate 
rotamer of one residue is in equations (13a,b,c), as follows. 
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( ^vdw 
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atoms / of side chain 
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(13a) 
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20 The molecular mechanics terms for a given rotamer may derive from averaging over 

all the sub-rotamcrs of that rotamer. 

The portion of the molecular mechanics energy measured in the reference structure i.e., 
VDW, HB and ELE is dependent only on the amino acid type, not its geometry, and thus is the 
25 same for each rotamer of a given residue. The role of the reference term is to help determine 
which sequences might be poor quality and not to distinguish between rotamer combinations 
of a particular sequence. The values of the reference energies in equations (13a, b,c) that are 
used in the molecular mechanics term are derived from calculations done with Perta over a 
large sample of main chain structures and sequences, whose results are averaged. 



30 



35 



The second term in equation (12), the main chain entropy cost, is also completely 
independent of the rotamer configuration and thus is an aid to distinguishing between 
sequences only. This term is computed just once for each sequence and is conveniently 
expressed as part of the intrinsic energy. The expression for the contribution to the 
main chain entropy from a given amino acid residue is in equation (14). 
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residue 

»t- type 



(14) 



The third term in equation (12), describes the side chain rotamer vibration entropy cost 
10 is measured with respect to a set of tabulated references. The weights of the sub-rotamere are 
obtained from the partition function. The contribution made by one rotamer to the side chain 
vibrational entropy component of the intrinsic energy* is expressed in equation (15). 
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(15) 
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In a preferred embodiment the vibrational contribution to the reference structure is 
20 calculated from a uniform distribution. In equation (16), N is the number of sub-rotamers, 
which can be adjusted by the user. 
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(16) 



In another embodiment, the reference term in equation (1 5) can be computed from a number 
3Q of calculations of Perla on a database of peptide fragments or, using a Gaussian Distribution. 

The fourth term in equation (12) which measures the solvation energy has been obtained 
by cutting the surface areas into intrinsic and pairwise parts. The contribution that one residue 
makes to the solvation component of the intrinsic energy is given in equation (17). 
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The solvation term is expressed from the relative buried surface area rather than the 
exposed surface area (thus the invention provides a subtraction in the sense of reference-target 
and not target-reference). For this reason, a different set of solvation parameters o f is used 
from those used to compute the solvation contribution to the scoring function energy of 
equation (1), as described later. 

Finally, the fifth term in equation (12) is a statistical contribution, which should consist 
of tabulated values, propensities or probabilities given by the user in a readable format. The 
contribution to this term in the intrinsic energy made by one residue is given in equation (18). 



~ ^Intrinsic^^stat ^ ^residue 
type 



(18) 
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5.4J. THE PAIRWISE TERM 



Finally, the pairwise energy term represents the interaction energy of a pair of candidate 
2Q retainers and is therefore summed over all pairs of candidate retainers. The previous comments 
about the reference states and temperatures are applicable here. 



The pairwise term comprises 4 contributions, equation (19): 



Pairwise = /r 




(19) 
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Similarly to the Intrinsic term, the molecular mechanics contribution to the pairwise term is 
split up into weighted contributions from van der Waals. electrostatic and hydrogen bonding 
energies. The summations run over pairs of atoms on different residue side chains, equations 
(20a,b,c). 
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(20c) 
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30 The molecular mechanics contribution to the pairwise energy also contains reference structure 
terms, i,e., VDW. ELE and HB 9 which depend only on amino acid types and are obtained by 
averaging over a number of different conformations obtained from a database. 

The second term in equation (19), for the rotamer vibration entropy, is formulated to 
^ measure the change of entropy due to the interaction of the two side chain ro tamers t ak i n g as 
a reference the vibration entropy of each rotamer substituted separately in the target structure, 
as shown in equation (21). 
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(21) 



20 



This entropy term is scaled by a factor X to avoid an overestimatian when summing over all 
25 pairs of interacting rotamers (see section 5.8 for details). 

The third term in equation (19), for the difference between accessible surface areas is 
formulated to measure the area buried between the two side chains. This solvation term, 
equation (22), is now also scaled by a factor X s to avoid an overestimation of the buried surface 
areas (likely to be counted several times when summing over all pairs of interacting rotamers). 
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The final term in equation (19) is. as previously, introduced to bias against improbable 
sequences of residues, though comprises a term to represent the prob ability of a pair f residues 
being present in the structure. The contribution for each pair of residues is expressed in 
equation (23), 

5 

~ >^P?iwise^^tat I* 1 ^rSdue ( 23 ) 
pair 
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5.5. THE MOLECULAR MECHANICS POTENTIAL 

In the method of the present invention, a protein is represented as an ensemble of atoms 
with discrete masses and partial charges and, therefore, classical mechanics equations are 
applied to estimate the potential energy of the system. 

The standard molecular mechanics function (or "force field") is a sum of terms that are 
related to bonded or nonbonded interactions and that depend on the atomic configuration, 
which is described by the coordinate vectors, r ( - (for an overview, see van Gunsteren & 
Berendsen, 1990, Angew. Chem. Int., Ed, Engl 29:992-1023), in equations (24a,b,c,d): 



20 



25 



4l- -'J- - *o) 2 *l^jK^ (O - Oof (24a) 
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Molecular mechanics is a well-established sphere of research and several widely used 
implementations exist: for example, AMBER, CHARMM, ECEPP2, MM2, CVFF, all of which 
are commercially or freely available. The operation of Perla is not dependent on the specific 
force-field which is used, because Perla is only concerned with ranking sequences. In an 
alternate embodiment, the user is free to adjust values of molecular mechanics parameters, 
where, for particular atom types, the parameters from a published force field are, in some way 
inadequate. 

10 



5,5.1, THE BQNP STRETCH ANP BOND-ANGLE TERMS 

The first three terms of the molecular mechanics force field correspond to bonded 
interactions. The first represents the elongation of covalent bonds between two atoms (bond 
stretching). It has a harmonic form, where b is the effective bond length, b 0 is the ideal length 
(energy minimum), and is the force constant that is characteristic of the actual type of 
covalent bond. The second term similarly describes the deformation of the angle <f> formed by 
three covalently bonded atoms (bond-angle bending). The third accounts for the rotation 
around bonds, or dihedral angles <p, according to a periodic potential with phase 6. 

The description of side chain conformations as a set of rotamers consists of setting the 
corresponding dihedral angles at values corresponding to energy minima. In addition, the 
covalent bond lengths and angles are set to their ideal values and are invariant. Therefore, in 
a preferred embodiment, the related energy terms are neglected and the methods of the present 

25 invention only consider the remaining three terms: van der Waals, hydrogen bonding and 
electrostatic interactions. These noncovalent forces that maintain protein three-dimensional 
structures, are the most important for a valid representation of protein structure, and are 
described in mathematical detail below. In a preferred embodiment, all parameters, e.£, , atomic 
charges and van der Waals energy parameters, ere taken from the ECEPP/2 potential (Momani 

30 et a/., 1975, /. Phys. Chetn. 79:2361-2381; Nemethy et aL % 1983, J, Phys. Chem. 87:1883- 
1 887). Although molecular mechanics force fields have been parameterized with goals other 
than protein design in mind, their application to protein design in the methods of the present 
invention can be justified because the overall scoring function contains other terms, such as 
entropic, solvation and statistical penalty terms, to account for respectively entropy, solvation 

35 and statistical bias, and because the contribution of the various terms can be adjusted by user- 
defined weighting factors. 
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5.5.2- THE VAN PER WAALS ENERGY 

Van der Waals interactions originate from a nonspecific attractive force that exists 
between atoms. That force is due to the transient asymmetry of the distribution of electronic 
charge around an atom, which induces a similar asymmetry in the distribution of electronic 
charge around neighboring atoms. The attraction increases as the distance between atoms 
decreases, until it is at a maximum when the two atoms, i and j\ are separated by a distance 
which is about 0.3-0.5 A larger than the van der Waals contact distance, (the closest contact 
distance between the two atoms that is observed in ciystal structures). The overlapping of the 
electron clouds of atoms i and j creates strong dominant repulsions at shorter distances (Fig. 

The van der Waals interaction energy between two atoms i and j can be described by 
a standard 6-12 Lcnnard- Jones potential, and the total van der Waals interaction energy term 
is the sum of the interaction energies between all nonbonded atom pairs: 



15 



■Evdw " Yd nonbonded// 



Aij By 1 



(25) 



where r & is the distance separating atoms i and j, and and are related to the type and 
20 environment of the interacting atoms. Thus, the energy term consists of a repulsive part that 
decays with r 12 , and an attractive part that varies inversely with r 6 . 

Van der Waals energies for pairs of atoms are on the order of the average thermal 
energy of molecules at room temperature (- 0.6 kcal/mol) and diminish rapidly even for a small 

25 increase of interatomic distances. Thus, the van der Waals interactionbecomes significant only 
when many interacting pairs form simultaneously, such as in the folded state of a protein. Most 
importantly, van der Waals energies are critical probes of the packing quality within the 
three-dimensional fold. For any sequence to fit a given fold, steric compatibility is required 
and no positive van der Waals energies should be tolerated. Cavities, which produce van der 

30 Waals energies near zero, should also be avoided, especially if they are small enough to exclude 
solvent water molecules. 

In the reference (denatured) state, atomic contacts within the polypeptide are less 
common and intra-molecular van der Waals interactions are not as significant as in the folded 
state, since van der Waals contacts with the solvent partly compensate for the loss of 
interaction. Therefore, most existing computer-based methods for designing proteins neglect 
the van der Waals contribution of the denatured state. 

34 
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However, consideration of the van der Waals energy of the reference state of a protein 
in the method of the present invention results in the removal of scaling artefacts in the van der 
Waals energy term. Potential energy functions that sum over all atoms of a system scale with 
the number of interacting atoms, resulting in energy terms that contain large numbers. These 
5 scaling artefacts that should be avoided. In the method of the present invention, scaling 
artefacts are avoided when comparing two sequences with different amino acid compositions 
by referencing each sequence to a denatured conformation. In the preferred embodiment of the 
present invention, van der Waals reference energies for each amino acid are utilized- These 
may be obtained in several different ways. In one approach, energy terms were calculated for 
10 each of the twenty amino acids in an extended five-residue peptide with alanine residues 
flanking the residue of interest. In another approach, energy terms are calculated for each of 
the amino acids, as found in a population of protein fragments similar to the population of 
unfolded structures. The reference energies scale with the number of atoms in each amino acid 
and compensate for the larger van der Waals contribution of larger residues in folded proteins. 
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5.5.3. THE ELECT ROSTATIC F.NFttfiV 

The interaction energy between electrostatic charges is fundamental and is described 
as the sum over all nonbonded atoms / and j, as follows, equation (26) 



^ec-Icharges/,^^ (26) 

25 

wherein q, and ty are the numbers of charges on atoms i and J, respectively, r tf is the distance 
between the two atoms, e is the charge of an electron, and e 0 and e, are the permittivity of the 
vacuum and the medium relative dielectric constant, respectively. 

3q In a vacuum, the electrostatic potential of an atomic charge in the field of another is the 

product of the two atomic charges divided by the distance that separates them (from Coulomb's 
Law). For two charges of opposite sign, the energy decreases as the atoms approach each 
other; the energy increases with decreasing distance if the charges have the same sign. The 
interaction is strong, e.g., up to 100 kcal/mol, and is effective over large distances. 

In media other than vacuum, the strength of the interaction is significantly reduced to 
less than several kcal/mol by the relative dielectric constant c r In one embodiment, the e, of 
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water has a value of about 80; the interior of a protein, which is mainly packed with carbon 
atoms, is less polar and has a lower dielectric constant, usually 4.0. 

In the method of the present invention, two dielectric constants (one for water or bulk 
solvent and one for the interior of the protein) are used for each mutated residue, according to 
the degree of burial of the side chain at the related position in the target protein structure. 
Every side chain atom is determined to be "non-core", i.e., exposed or "core", i.e., buried 
according to some geometric criterion. In one embodiment, this criterion is derived from the 
relative proportion of the solvent accessible surface area of the side chain of the input structure 
that is assigned to the atom. In a preferred embodiment, the geometric description takes into 
account the distance from Ca to the nearest solvent molecule on a solvent accessible surface 
constructed with an 8 A probe, taken along the Ca - Cp vector, as well as the shortest distance 
from the Ca atom to any solvent molecule on that surface. For each pair of atoms for which 
an electrostatic interaction is being calculated, the dielectric constant used will be the solvent 
1$ value if both atoms are exposed, the protein interior value if both atoms are buried. When the 
interaction is between a completely buried atom and a completely exposed atom, the average 
of both dielectric constants is used. 
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In addition, the electrostatic energy term is modified to lessen the importance of the 
interaction at long distance. In one embodiment, shown in Equation 27, dielectric constants 
are scaled linearly with the separation distance between atoms i and j: 

z r = zr i} (27) 

In the preferred embodiment of the present invention, the pH is considered to be neutral, and 
the parameters used are for fully charged versions of acidic (aspaitic acid, pK^.S; glutamic 
acid, pK«4.5; histidine, pK=6) and basic (lysine, pK=ll; arginine, pK-12) amino acids. 
Therefore, the entire electrostatic energy term is scaled by an exponential factor to account for 
the screening of charges by salts and counterfoils, as shown in Equation 28: 
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35 The rate of exponential damping is controlled by a decay constant, r t , whose units are those 
of distance. The decay constant, r 3 in equation (28) can be calculated according to methods 
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described in Lacroix E., Viguera A.R. and Serrano L. 1998 J. Mol. Bi 1. 284:173-191. An 
expression for r, is given in Equation (29): 



8fte 2 N A 7 
1000AT 



(29) 



In equation (29), N A is Avogadro's number, / is the ionic strength of the solution and k is 
10 Boltanann's constant. 

When the electrostatic interaction energy is computed according to equations (26) 
and (27) or (28) the preferred value of e r is between 8.0 and 32.0. Fig. 7 illustrates the 
electrostatic interaction energy between two unit charges of equal sign as a function of 
interatomic separation calculated using either Equations (26) and (27) or equation (28). 

In the denatured state, electrostatic energy terms contribute less to the potential 
energy due to the overall increase in interatomic distances caused by extension of the 
protein chain, and more importantly, to higher solvation and charge screening. In contrast, 
in the folded protein, amino acids separated by only a few residues rarely undergo a 

10 conformational change that gives rise to a significant change in the distances separating 
their atomic charges. Therefore, if only the native state of the protein is considered, the 
electrostatic term is over-estimated. In a preferred embodiment of the present invention, 
values are tabulated to represent all possible electrostatic interactions in the denatured state 
as a function of the sequence separation. In one embodiment of the present invention, the 

! * electrostatic energy term of the reference state is zero. 



5.5.4. HYDROGEN BONDING ENERGY 

A hydrogen bond is formed when two electronegative atoms, a donor and an 
acceptor, compete for the same hydrogen atom. As a result, the distance between the 
hydrogen atom of the hydrogen bond donor and the hydrogen bond acceptor is shorter than 
the van der Waals contact distance, although it is larger than the length of a covalent bond. 
The interaction is partly covalent and partly electrostatic in nature and can have an energy 
ofupt 7kcal/mol. Hydrogen bonds are highly directional and occur predominantly with 
the donor, hydrogen, and acceptor in a collinear orientation. Therefore, the potential energy 
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function of the preferred embodiment of the present invention considers hydrogen bonding 
only if the geometrical conditions are satisfied, ie. 9 if the distance between the hydrogen 
and the acceptor atom is between 1.7 A and 2.5 A, and the angle made by the donor, 
hydrogen and acceptor is greater than 100° (Fig. 4). If these conditions are met, a hydrogen 
bonding term (Equation 30) replaces the van der Waals teim corresponding to the 
interaction between the hydrogen and acceptor atoms. 
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(30) 



The preferred embodiment of the present invention does not take into account the 
possibility that there is intra-molecular hydrogen bonding in the denatured state, since the 
geometrical conditions are only fulfilled if elements of structure, e.g., turns or a-helices, 
15 form locally. The formation of hydrogen bonds between atoms in the denatured protein and 
water is included empirically in an accessible surface-area-dependent solvation potential 
described below in Section 5.6. In essence, the residues in a denatured protein are modeled 
by ensembles of representative fragments taken from protein structures in the PDB. 



20 

5.6. THE EMPIRICAL POTENTIAL 

The method of the present invention, as implemented in the preferred embodiment, 
Perla, also evaluate changes in entropy and solvation of the protein chain by means of 

25 empirical models constructed to account for properties that cannot be broken down into a 
set of well characterized physical forces. It is customarily difficult to accurately model 
entropy and solvation, because a practical and accurate representation of the ensemble of 
unfolded protein structures would have to be developed This would necessitate the 
handling of an enormous number of either water molecules or chain configurations within a 

30 practical amount of computing time, as well as the development of an accurate set of energy 
parameters to describe these unfolded states. Instead Perla adopts pragmatic levels of 
approximation. 
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5.6.1. SOLVATION 



Proteins function in aqueous media, which are poor solvents for apolar molecules 
because apolar molecules cannot participate in hydrogen bonding with liquid water. To 

5 satisfy their hydrogen bonding requirement, water molecules that surround a hydrophobic 
molecule order themselves by hydrogen bonding with each other, and consequently, lose 
many degrees of freedom. The reduction of exposed hydrophobic surfaces through protein 
folding leads to a release of ordered layers of water molecules, and consequently, the 
entropy of the solvent increases. This increase in the entropy of the system is the basis for 

1Q the hydrophobic effect, which leads to proteins adopting compact shapes. In terms of 

protein design, the essential property of water is the partitioning of polar and apolar residues 
between the protein surface and interior, or core. As a result of the hydrophobic effect, 
apolar residues arc preferentially, but not always, buried in the protein interior, where the 
aqueous solvent is excluded. Conversely, polar residues may occasionally be buried but are 

15 preponderantly found on the protein surface; charged residues are rarely buried. 

Due to the large number of water molecules in the layers surrounding the protein 
surface, water cannot be explicitly modeled in order to consider the effect of solvent on 
sequence preferences. Eisenberg and McLachlan showed that the free energy of interaction 
of a protein with water could be represented by the sum of the interaction energies of each 
atom of the protein with solvent. They further proposed that the interaction strength is 
proportional to the accessible surface area, (ASA,), of each atom. 

In a preferred embodiment of the present invention, water is modeled implicitly (in 
bulk, rather than as discrete molecules) and the solvation potential energy term is calculated 
25 from the difference in accessible surface area of each atom / in the folded and denatured 
protein (A ASA,) and from empirically determined solvation parameters for each atom (a), 
as shown in Equation 31 : 



30 



Accessible surface areas depend only on the atomic configuration of the protein and 
are calculated using the method of Lee and Richards (1971, J. Mol Biol 55:379-400) and 
the numerical surface calculation (NSC) routine of Eisenhaber et al. (1995. / Comp. Chem. 
16:273-284). Many other suitable surface area calculation algorithms are. however, known 
to one skilled in the art and could be utilized with the methods of the present invention. A 
water molecule with radius of 1.4 A is the "probe" that is rolled along the van der Waals 
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surface of the protein atoms in order to calculate the accessible surface. The atomic radii 
and solvation parameters used to calculate the accessible surface areas of proteins in a 
preferred embodiment of the present invention are taken from Eisenberg and McLachlan 
(1986, Nature 319:199-203), as described in section 5.11. 

However, correct accessible surface area measurements can only be made in the 
context of a full structure, and not before the optimal combination of side chain rotamers is 
found for the evaluated sequence. Therefore, in order to evaluate solvation contributions to 
the intrinsic and pairwise energies, an alternative approach is adopted. 

10 It has been found that, when evaluating the intrinsic and pairwise solvation energies 

for different sequence variants and conformations during optimization, changes in ASA 
values for pairs of residues upon folding leads to an overestimation of the area of buried 
surfaces. Nevertheless, it is important to include a solvation parameter during optimization 
routines. Thus, in a preferred embodiment of the present invention, polar and apolar buried 

15 

surfaces evaluated in a pairwise manner are scaled down in the manner proposed by Street 
and Mayo (1998, Folding <fe Design 3:253-258). The surface area of residue i buried by the 
template, given by BSA h is evaluated as the difference between the accessible surface area 
of the same residue placed at the center of a five residue peptide and the accessible surface 
area of the residue in the target protein conformation, is shown in equation (32): 

20 

target 

BSJk = ASAf-' 6pddfS - ASjC"™ (32) 

The surface area buried between residues i and j is evaluated as the difference 
^ between the exposed surface area of each residue separately placed in the target 

conformation and the exposed surface area of the pair of residues placed together in the 
target protein conformation. This is shown as follows: 



30 
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target 
structure 



target 
.structure 



target i 
structure 



+ ASAj - ASAiwd J 



(33) 
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To compensate for the overestimation of total buried surface area, the ASA terms are scaled 
with a factor, X„ that depends on the location of residues i and/ Equation (30) differs from 
that given by Street and Mayo because the factor X 5 multiplies all the ASA terms. In one 
embodiment, X t is taken to be 0.40 for core residues, 0.75 for non-core residues, and 0.60 
5 for a pair that consists of one core residue and one non-core residue. In a preferred 
embodiment, X, can be related to an alternative, pre-calculatcd parameter, A: 



^Contacts , xrcontacts 

-"first rotamer 1 v second rotamer ( 

**s = A ^contacts ^contacts ^ ' 

1 0 iV first rotamer iV second rotamer 



In a preferred embodiment, A « 0,5, though others values could be obtained by a 
fitting exercise. The solvation energy is obtained by summing equations 32 and 33 over all 
side chains and by multiplying the accessible surface areas by the solvation parameters 
0.100 for polar buried surfaces and -0.026 for nonpolar buried surfaces. 



S.6.2. ENTROPY OF THE MAIN CHAIN 

The entropy change upon folding, another major component of protein stability, is 
calculated in a preferred embodiment of the present invention using a statistical approach- 
Although entropy is a unified physical concept, it is practical to divide the entropy change 
into parts related to either the main chain or to the side chains (see Section 5.9, below). The 
main chain entropy term is expressed as the cost to fix the backbone conformation into the 
ensemble of <J> and \Jr dihedral angles of the target structure. These costs depend on the 
nature of the amino acid located at each <J>-iJf pair, and were predetermined for use in the 
preferred embodiment of the invention to reflect the secondary structure propensities of the 
twenty amino acids (Mufioz & Serrano, 1994, Prot eins 20(4):301-311), as described below. 
3Q The entropy of the main chain can therefore also be thought of as a secondary structure 
propensity term. 

A set of 527 protein structures that share less than 35% sequence homology 
(PDBSELECT; Hobohm & Sander, 1994, Protein Set 3:522-524; Hobohm et al 9 1992, 
Protein Set. 1 :409-417) was used to obtain all main chain dihedral angles. The numbers of 
occurrences of each amin acid in regions of the Ramachandran ((|>-t|0 ptot sampled at fixed 
intervals, d Q were determined. In a preferred embodiment, d Q is taken to be twenty degrees. 
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The tendency for amino acid X to populate a particular region of Ramachandran space, e.g., 
fyrtyi, is the ratio of the number of hits in the interval considered (N^ 4 ) and the total 
number fhits for amino acid X(N iH ^): 

x s I±fL. (35) 

For such a partitioning of dihedral angle space, it is also necessary to quantify the distance, 
d 9 between pairs of points (each of which represents a residue conformation): 



10 



^{✓20°Jf20° = Vfe " hQ°X20 o f + (Vt ~ ^20°Jf20°) ( 36 ) 



15 



For 20° X 20° regions which are more than a threshold distance (in angle space), d t from 
the residue <j)j - % dihedral angles, the occurrences are modified with a weight given by an 
exponential decay function of the separation distance, (to the inventors* knowledge, such an 
20 approach has not been used before in computer-based protein design methods): 



^20**20* = ex P("^^20°^20° ,d o ) (37) 

25 

This modification allows a smooth transition of the energy function over the main chain 
dihedral angle space (instead of the abrupt changes that occurred previously when crossing 
the 20° x 20° boundaries). Thus, empirical observation is used to calculate the likelihood 
20 that a particular amino acid will have main chain dihedral angles 4>, and i|r, in any given 
protein structure. 

In the preferred embodiment, the database used is large enough to be considered as a 
system under thermodynamic equilibrium in which pseudo-energies or costs to displace the 
equilibrium toward a particular state can be calculated from the natural logarithm of the 
35 ratio in Equation 32. Entropy costs to fix the main chain dihedral angles f amino acid X in 
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a particular region of the Ramachandran plot, e.g., 4> r iJ/„ were therefore calculated as shown 
in equation (38) for use in a preferred embodiment of the present invention: 



10 



residue 

Eat type 

subspaccs ^^20^20° 

« close to residue , 7CA 

' KI phys ln residue W 



5.7. CONSIDERATION OF THE DENATURED STATE OF PROTEINS 

1 5 An important consideration when modeling proteins is that of a reference state or 

configuration. In practice, proteins in solutions are dynamic ensembles of structures: the 
folded (•'native*') structures are in equilibrium with unfolded "denatured" configurations. 
The former are compact, with residues in close contact with non-neighboring residues due 
to the intricacies of the backbone configuration, whereas the latter are open, more chain- 

20 like. For the purposes of the present invention, the essential difference between these two 
extremes is that individual residues (particularly their side chains) participate in many more 
pairwise interactions amongst themselves in the folded state than in the unfolded. It is 
desired to quantify this difference at the level of individual residues, a result which is 
achievable as described below. A further reason to use a reference state configuration is s 

25 that the calculated solution store remains approximately proportional to the stability of the 
solution structure. 

Whereas native protein structures are available (in the PDB) denatured structures 
must be obtained via simulation. In a preferred embodiment, a set of non-homologous 

30 proteins (obtained from the WHATIF database) is used to extract all protein fragments that 
are at least 4 and at most 20-residues long. These peptide segments may be clustered into 
groups according to length and structural homology, using a combination of main chain 
dihedral angle comparisons and internal (i.e.. Ca-Ca) distance comparisons. There are 
many clustering algorithms which may be used for this purpose, for example, Ward's, 

35 Jarvis-Patrick and assorted hierarchical methods. For each cluster, a single representative is 
selected (for example, from the geometrical center of the cluster). The ensemble of 
representatives is used as a set of main chain templates to reconstruct sequences of the type 
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(Ala) M -X -(AIa)„ and (Ala) r X-(AIa) m -Y-(Ala) fl where X and Y are any of the 20 natural 
amino acids (and the subscripts, /, m, n, represent segment lengths) and Ala represents the 
amino acid, alanine. (These sequences represent the amino acid residues of interest in an 
ordinary environment: alanine is the amino acid with the smallest (shortest) carbon 

5 containing side chain. It therefore contains no polar or bulky groups which might cause 
folding or twisting of the sequence but, unlike glycine (whose side chain is hydrogen), has 
sufficient bulk to prevent collapse.) Perla itself can be used to determine a solution score 
for each sequence, i.e., each peptide fragment. It is then possible to compute an average 
solution score that corresponds to the output of a partition function measured over the 

10 ensemble of fragments. Perla does so for each separate energy term, and then provides sets 
of values to be used as reference values for the random coil, i.e. 9 the denatured state. The 
references for the intrinsic parts of the van der Waals, hydrogen bonding and electrostatic 
energy terms, and the side chain rotamer entropies, are measured with sequences of the type 
(Ala) m -X-(Ala) { , while the references for the paiiwise parts of the van der Waals, hydrogen 

j 5 bonding and electrostatic energy terms, are measured with sequences of the type (Ala)p X- 
(Ala) m -Y- (Ala) ft . In a preferred embodiment of the present invention, when more than one 
CPU is available, reference energies may be computed in parallel with the energies of the 
folder structures. 

For application to other categories of macromolecules, it may be appropriate to 

20 

consider an alternative form of reference state. For example, when attempting to find a set 
of nucleotides that improves the interaction between a nucleic acid and a protein or other 
nucleic acid sequence, the reference may be the isolated nucleic acid in question, whereas 
the scoring function will quantify the extent of the interaction. 



5.8. CONSTRUCTION OF THE ROTAMER LIBRARY USED IN PERLA 

Crystallographically determined protein structures show that the side chain dihedral 
angles are not distributed uniformly through 360° (Janin & Wqdak. 1978, J. Mol Biol 
125:357-386; Ponder & Richards, 1987, /. Mol Biol 193:775-791), For example, the 
torsion angles around two bonded sp 3 carbons generally cluster into 'gauche + 1 (+60*), 
'gauche - (-60°), and 'trans' (180°) conformations (Fig. 6). In a preferred embodiment of 
the present invention, we wish to construct rotamer libraries in which all significant 
confonners are represented. 
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By way of example, the set of 527 protein structures that share less than 35% 
sequence homology (see Section 5.5; PDBSELECT; Hobohm & Sander, 1994, Protein ScL 
3:522-524; Hobohm et al, 1992, Protein Sci. 1:409-417) was used to obtain all sets of side 
chain dihedral angles (x„ %n X*> 7U)- For each amino acid other than glycine and alanine, 
5 the distribution u(x) of dihedral angles determined from the Protein Data Bank structures 
was fitted to a combination of normal distributions represented by a sum of Gaussians, 
equation (39) as follows: 



10 



( J ^ 



exp 



(-(*-*) 



2\ 



Zcxf 



(39) 



15 



The number of Gaussian terms, N, was modified until no further reduction of the square of 
the difference between observed and calculated distributions was obtained, A constant 
term, k, was added to fit the distribution of side chain dihedral angles with poorly marked 
preferences, e.g., x 2 of Asn and Asp, or to represent the noise. The outputs of the fitting 
procedure are the centers (ji r ) of the N Gaussian peaks and their standard deviations, o.. In 
addition to the expected well-defined peaks of standard conformers, separate distributions 
of lower amplitudes were used to fit the data. In most cases, as illustrated for valine (light 
grey peaks in the top panel of Fig. 12), the additional Gaussian curves represent variations 
around different dihedral angle values that can be adopted by a particular amino acid in a 
significant number of instances. 

Side chain rotamers with all combinations of the peak centers, \i t (except for "ghost" 
peaks) were constructed using ideal values for the covalent bonds and angles (Mazur & 
Abagyan, 1989, J. Biomol Struct. Dyn. 6(4):8l5-832; Momani etal, 1975, /, Phys. Chern* 
79:2361-2381; Nemethy et al, 1983, J. Phys. Chem. 87:1883-1887). For dihedral angles 
that do not have a normal distribution, e.g., % 2 of asparaginc and aspartic acid, and Xj of 
glutamine and glutamic acid, sets of values were chosen to sample the range of observed 
values. The constructed side chains were inspected visually and conformations with steric 
clashes were removed. The remaining side chain conformati ns form the custom-made 
rotamer library of 419 rotamers employed by the preferred embodiment of the present 
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invention, (Table 1). Side chains built by the preferred embodiment of the present 
invention are taken from this library. The advantage of this approach for the construction of 
a rotamer library is that it does not use stereochemical rules to generate the retainers, thus 
allowing for the addition to the library of less abundant but relevant retainers. Furthermore, 
5 the fitted normal distributions define the margins within which the rotamers can oscillate 
during evaluation of sequences. 

The generation of dihedral angles by means of equation (39) does not include the 
fact that consecutive angles have correlated distributions. The correlation, due to the 
topological structures of certain amino acids, can be so strong that a particular value of X\ is 
10 only possible if x 2 itself adopts a defined value (e<g, , the -95, 36 conformer of leucine in 
Fig. 12). Furthermore, the fitting procedure used to generate the dihedral angles for the 
rotamers used in the library cannot detect rare peaks. 

More preferably, side chain conformations with correlated dihedral angles or rare 
1 S dihedral angles can be included in the rotamer library by repeating the analysis above while 
considering the distribution of all dihedral angles of an amino acid simultaneously, since 
rare peaks are more resolved in a multidimensional representation. The rotamer library 
employed by the methods of the present invention contains only a few identified cases of 
rare or correlated dihedral angle values, e.g., the -175, 150 and -145. -150 leucine 
20 conformers, (Fig. 12). 



Table 1 : Summarized Description of Rotamer Library 



25 



Amino Acid 


Dihedral 


Number of 


Amino Acid 


Dihedral 


Number of 




Angles 


Rotamers 




Angles 


Rotamers 


Ala 




1 


Met 


Xit X2> Xi 


27 


Cys 


Xl> Xz 


18 


Asm 


Xi.Xz 


18 


Asp 


Xi» Xz 


9 


Pro 


Xi. X2> Xy 


3 


Glu 


X\> Xz> X) 


27 


Gin 


Xi» Xz» X3 


54 


Phe 


Xi> Xi 


9 


Aig 


Xlr Xz» Xi> Xa 


75 


Gly 




1 


Ser 


Xi» Xz 


18 


His 


Xi» Xz 


12 


Thr 


Xi» Xz 


18 


lie 


Xi.X 2 


9 


Val 


Xi 


3 



46 



BNSOOCIO: <WO 0) l88tOA2_L> 



WO 01/16 



PCT/EP00/08504 



Lvs Z1.X1.fc.X4 77 Tip X |.fc 12 

Uu ' X..fc 10 Tyr X|tXa 18 

TOTAL: 419 

5 

The use of retainers enables computation of side-chain vibrational contributions to 
the entropy. The calculation of side chain entropy terms is described subsequently in the 
discussion of the mean field approximation. 

I o The intrinsic vibrational entropy term represents the change from a uniform 

distribution to the distribution obtained from the partition function described as the equation 
of Section 5.8.2, and the pairwise vibration entropy term is the change from the initial and 
main chain-dependent distributions of sub-rotamers of the two side chain retainers to the 
distribution of sub-retainer pairs due to their interaction with each other. Scaling of the side 

! 5 chain vibrational entropy contribution to the pairwise energy by a factor X is necessary to 
avoid the overestiraation of the entropy change when summing over all pairs of interacting 
retainers. This arises for a similar reason to the over-estimation of the solvation energy. In 
the preferred embodiment of the present invention, A. is given by equation (40): 

20 

kf-contacis , ^contacts 
first rotamer "** second retainer 



i = A 



xr contacts ^contacts (40) 
iV first rotamer iV second rotamer 



25 

Here, A is set to be 0.5. 



3Q s -5- ENERGY MINIMIZATION AND ELIMINATION 

OF INCOMPATIBL E AMINO Arm CONFORMERS 

There are two facets of the energy minimization: minimization of the interaction 
between side chain and template; and niinimization of the pairwise interactions between 
side chains. In either case, there are at least two possible methods of minimization. 
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5.9.1. METHODS OF ENERGY MINIMIZATION 



Sidd chain conformations taken from the rotamer library of the preferred 
embodiment of the present invention are idealized. Consequently, in the preferred 



embodiment of the present invention, their interactions with the protein template and other 
* side chains are optimized and flexibility is introduced to relax strain inherent in the fixed 
geometries provided by the rotamer library. 

Energy minimization is carried out in dihedral angle space using the non-bonding 
terms of the molecular mechanics force field (van der Waals, electrostatic, and hydrogen 
10 bonding). 

In one embodiment of the present invention, the energy minimization method may 
be so-called "Steepest descent"; in another, it may be taken from a class of methods known 
as "quasi-Newtonian* 9 . The theory behind these methods is accessible to one skilled in the 

^ j art (and can be found in Numerical Recipes in C - The Art of Scientific Computing by WH 
Press, 2 nd edition, section 10.6, SA Teukolsky, WT Vetterling and BP Flannery), but in 
general terms, the interaction energy and its gradient with respect to displacements in 
dihedral angle space are utilized Minima on the energy surface are located iteratively 
through use of the gradient to search downhill from a given point. The quasi-Newton 

2Q methods attempt to gather information about the curvature of the energy surface. Methods 
in this category include "BFGS" and "conjugate gradient", the distinctions between them 
arising in how each decides to approximate the Hessian (matrix of second derivatives of the 
energy). In practice, the method of conjugate-gradient has been found to be effective, 
provided that certain precautionary measures are taken in order to avoid large rotations that 

25 would in fact transform one side chain rotamer into another one (this would deteriorate all 
subsequent partition functions, i\e, that calculated to eliminate the rotamers that have the 
lowest probabilities according to the interaction energy with the main chain). First, the 
rotation step size has to be small (fractions of a degree to a few degrees); thus the factors 
that multiply the gradients are small, and the gradients themselves are capped at some 

30 energy maximum. Second, the energy function is modified to contain a rotation penalty 
function to directly limit the minimization sampling to conformations close to the initial 
rotamer structure; it has the form in equation (41). 



35 




(41) 



minimized^ 
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The "force constants" k z are expressed as functions of the standard deviations measured 
during the construction of the rotamer library by fitting of the frequency distributions 
observed in the protein database as given by the o- values in equation (39) for example. 
Third, if a large rotation (more than a couple of standard deviations) is done despite the 
5 penalty, the rotamer is simply placed back in its initial conformation, discarding the result 
of the minimization. 

In a preferred embodiment of the present invention, minimization is carried out by 
exhaustive sampling of dihedral angles of rotamers close to the ideal conformation. This 
method is not only simpler, but is superior to conjugate-gradient and similar methods of 

10 

optimization. (As an example of why this is so, the formation of a hydrogen bond may 
require an energetically unfavorable rotation of a side chain in order for the correct 
geometry to be achieved, which would only be discovered using a method that samples 
rotamer conformations close to the ideal conformation without energy minimization,) The 
conformations which are obtained through systematic rotations around the dihedral angles 
X, and Xi are referred to as "sub-rotamers" In a preferred embodiment of the present 
invention, the step size and number of steps are precompiled for each of the twenty 
naturally occurring amino acids and are determined to cover rotations smaller than two 
standard deviations away from the minima in dihedral angle space (derived during the 
creation of the rotamer library). Such a range is usually about 1 S degrees, or even smaller 

20 

than a single standard deviation if it is necessary to optimize the residue set. It is expected 
that a sequence which enables the packing of rotamers in their ideal conformation (that of 
the library) should be preferred to another sequence that would necessitate rotations of its 
side chain rotamers. Although there is an advantage in using many small steps in the 
thoroughness of coverage, the calculation time has to be considered. In the preferred 
embodiment, three steps of 5 degrees in each direction around the dihedral angles give 
good results, i.e., 7 steps in all for each angle. This leads to 49 sub-rotamers, and 49x49 
energy calculations to obtain a pairwise energy term. 

When considering the interaction of the side chain with the template, the final 
30 energy is the weighted average over all possible sub-rotamers. The partition function that 
defines the weight of state / as a part of a system with N states with energies E, is defined as 
follows: 



35 
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j(includingi) 



Accordingly, the contributions of sub-rotamers to the intrinsic energy during 
10 minimization are given by equation (42). The E t comprises molecular mechanics energy 
terms. 



15 



In a preferred embodiment of the present invention, pairs of side chains are 
minimized using systematic sampling of sub-rotamers. The outcome of using sub-rotamers 
to optimize pairs of side chains is the weighted average of all possible pairs using a formula 
such as (43): 



exp 



20 



w 0 = 



RTJ 



EZ ex P 

' J 



RT 



(43) 



Minimizing pairs of side chains individually (regardless of other side chains) is 
assumed to be correct as long as the actual conformations of the side chains depart only 
slightly from their starting point (to prevent any incompatibility in larger sets of side 
chains). 

30 

5.9.2. PROCESSING RESULTS OF MINIMIZATION 

In the method of the present invention, after the intrinsic energy computation, side 
chain confonners that do not interact favorably with the protein target conformation after 
3S minimization are rejected. In one embodiment of the present invention, retainers for which 
the intrinsic energy term is above a predetermined threshold are rejected. The use of an 
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absolute threshold, however, may not be ideal. Intrinsic energies vary in magnitude from 
site to site because they depend on the actual backbone structure. Poorly resolved structures 
tend to show many spots with local repulsions. Moreover, strains do exist in proteins where 
repulsions are common, e.g., in turns, where the tight reversal of the main chain in a turn is 

5 often accomplished by placing the <J>-\|r dihedral angles in less favored regions of the 
Ramachandran plot. Therefore, in a preferred embodiment of the present invention, the 
absolute energy threshold is placed high enough (about 50 kcal/mol) to accept enough side 
chain confonners at every position of the target structure, and a subsequent relative 
threshold is then applied to keep only the most qualified rotamers. The subsequent relative 

10 energy threshold is designed to exclude those rotamers whose partition coefficient weights 
are very small The absolute energy threshold is scaled according to the various user- 
defined weighting coefficients. When all w are 1.0, 50 kcal/mol" 1 is sufficient, but when all 
w=0.5, except for solvation, a value between 10-20 kcal/mol" ! , is appropriate. The 
subsequent relative threshold is determined by calculating for each amino acid a partition 

j 5 function oyer all rotamers using the template-side chain interaction terms, equation (42) 
(i.e., an intrinsic energy term, summed over the remaining rotamers) with the minimum 
threshold being fixed as a minimat probability of frequency, e.g. 9 0.001 ♦ 



20 



25 



30 



35 



The way in which minimization may be applied to other categories of 
macromolecules will depend upon the complexity of the building blocks and upon the 
overall structure of the macromolecule. In long extended systems such as nucleic acids, the 
relevance of the painvise energy term will be less important. Similarly, in systems with 
inflexible sidechains or with little conformational flexibility, the need for a thorough 
minimization protocol will be diminished. 



5.10. OPTIMIZATION ROUTINES 

For peptides longer than a few residues, an exhaustive sampling of every possible 
sequence and combination of rotamers is not practical. Hence, in the methods of the present 
invention, procedures are used which either decrease the size of the sequence space to be . 
covered or rapidly find the probabilities of having particular rotamers at each position of the 
modeled protein structure. These procedures are often termed "semi-exhaustive" 



5.10.1. DEAD-END ELIMINATION 
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In ttie method of the present invention, significant reductions of sequence space are 
obtained by discarding amino acids that cannot belong to the optimal sequence, which is 
that with the lowest potential energy, and the calculation time is shortened. To eliminate an 
amino acid, the modeling procedure has simply to exclude all its known side chain 
s conformations. Undesired retainers are detected by applying the dead-end criterion, which 
is the underlying principle of the dead-end elimination ("DEE") theorem. The theorem 
states that, for a given residue j. a particular rotamer, i„ is not compatible with the global 
minimum energy conformation ("GMEC") if, for the same residue i, an alternative rotamer 
i, exists for which the following inequality holds true (Desmet ei a/., 1992, Nature 356:539- 

10 542 ) : 



' j*i r 1 J** 

15 The minimum and maximum functions (min, and maxj cycle over all rotamers/, of 

residues j, searching for the rotamer which offers, respectively, the lowest and highest value 
of the interaction energy with residue i. The rotamers picked by the minimum or maximum 
function, j„ do not have to be identical. Thus, the left-hand side of Equation 44 evaluates 
the best possible interaction of rotamer i r with the side chains of all other modeled residues, 

20 while the right-hand side evaluates the worst possible interaction of an alternative rotamer i, 
for the same residue with all the other modeled residues. A side chain rotamer is "dead- 
ending" if its best interaction with the surroundings is less advantageous than that for 
another rotamer of the same side chain taken at its worst. Only one rotamer i, that satisfies 
the inequality has to be found to qualify as dead-ending. 

25 In another embodiment of the invention, a more powerful version of the elirnination 

criterion is utilized that is less restrictive and therefore more effective (Goldstein, 1994, 
Biophys. J. 66:1335-1340). It states that a ride chain rotamer i r is dead-ending if the energy 
of the model can be lowered by the choice of an alternative rotamer i p while keeping all 
other side chains fixed. This elimination criterion is described in Equation 45 for the same 

30 set of/,: 

^template _ ^template + £ ^(^ise _ ^pairwise j > Q ^ 

35 
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Both dead-end elimination criteria can be extended to pairs of rotamers. The energy 
contribution of a rotamer pair (ij t ) and its interaction with a third residue/rotamer k, ate 
given by Equations 46 and 47, respectively: 

^template -template -pairwise 



10 



E (f,M ■ *ftA) + E Um*i) <" > * *> < 47 > 



Then, Equation 44, in one embodiment of the present invention, can be written for paiis f 
^ rotamers, as follows: 



E 0, Jt ) + X "»faf .;,).*, ) > %.^> + ^ msx <( E «M* t ) (48) 

20 



and in a preferred embodiment of the present invention, Equation 43, can be rewritten for 
pairs of rotamers as follows: 

25 

%*) - *<wv) + £ "Mfyj.** - ) > 0 (49) 

30 

Dead-ending pairs do not lead to an elimination of a particular amino acid unless 
one of the participating rotamers is the only possible side chain conformer for the related 
residue position, in which case die other rotamer of the pair is not compatible with the 
GMEC and can be discarded. Lasters and co-workers (Lasters et al, 1995, Protein Eng. 
35 8:81S-822; Lasters and Desmet, 1993, Protein Eng. 6:717-722) showed that dead-ending 
pairs can be safely ignored in the minimum term of Equation 44 and the left-hand term of 
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the minimum function of Equation 45. Due to the exclusion of dead-ending pairs, the 
minimum functions might return higher values that strengthen the rotamer elimination 
criterion. 

In the preferred embodiment of the present invention, the dead-end elimination 
5 routine follows an iterative process as follows; (a) dead-ending rotamers are eliminated, 
repeating evaluations of Goldstein's criterion (Equation 45) until no more dead-ending 
rotamers are found; (b) dead-ending pairs are detected using the first elimination criterion 
(Equation 44; Desmet et ai, 1992, Nature 356:539-542) because it is estimated more 
quickly; and (c) new cycles of rotamer removal as in step (a) are carried out. This continues 
10 until no more dead-ending pairs can be found. The more effective but computationally 
expensive criterion for the detection of dead-ending pairs (Equation 48) is used when many 
rotamers have been eliminated and the whole cycle is restarted. At the end, if all rotamers 
of an amino acid at a particular site in the protein are dismissed, sequences containing this 
particular amino acid at that site are also abandoned. 

15 

In one embodiment of the present invention, the DEE routine can be used to 
determine an optimal set of rotamers for a given sequence. In a preferred embodiment, the 
routine is not used to limit the output to one optimal set of rotamers. since side chains, 
particularly those positioned on the solvent-exposed surface, are flexible and adopt different 
20 configurations. 



5.10.2. MEAN FIELD THEORY 

25 The foregoing will have provided the user with one or more acceptable rotamers for 

each residue position. This outcome represents a level of statistical uncertainty in Ihe 
situation being described. It is not adequate to subsequently model the system with merely 
one rotamer for each residue. 

It is desirable to obtain a contribution to the entropic term from all possible 

30 

conformations of side chains that remain after iterative dead-end elimination has eliminated 
all sequences that cannot belong to the global energy minimum, as described above in mis 
section. In a preferred embodiment of the present invention, mean field theory (MFT) is 
utilized to achieve this. MFT is an iterative technique which has found wide application in 
35 the physical sciences for describing systems of interacting particles which may adopt many 
different energy states. 
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All possible side chain conformations are considered using a mean field 
approximation that is designed to provide an estimate of the entropy of side chains in both 
the denatured and the modeled state, i.e., the template. The method attributes to each side 
chain conformation of all residues in the protein sequence that are not fixed a probability 
that depends on the average of all possible environments, weighted in turn by their 
respective probabilities of occurrence (Koehl & Delarue, 1994, J. Mol. Biol. 239:249-275; 
Koehl & Delarue, 1995, Nat. Struct. Biol. 2:163-170; Kochl & Delarue, 1996, Curr. Opin. 
Struct. Biol. 6:222-226). The probability of occurrence is related to the energetic 
favorability. The system is initialized by giving ah equal probability to each side chain 
rotamer so that every side chain conformation "feels" equally the presence of all rotamers at 
other residue positions. At this point, the field energy of each rotamer is the sum of all 
possible pairwise interaction energies normalized by the number of interactions, plus a 
contribution from the interaction with the protein template, as shown in equation (50): 

all residues all rotamers 
j * 



20 



25 



30 



wherein MF(i t ) is the mean field energy experienced by rotamer r of residue L Initially, the 
probabilities, w, for the rotamers of any residue are equal. At all times, in order for them to 
be interpretable as such, these probabilities sum up to 1. The application of MFT is to 
minimize the terra AfFby suitable adjustments of the weights. 

Having obtained the mean field energy perceived by each rotamer of a residue, 
proper weights can be estimated from a partition function that integrates all field energies, 
so that the rotamer that interacts best with its environments becomes more probable than 
competing rotamers equation (48): 

l-MFA 

w, r * ~~k — n^T (SI) 



all rotamers/ 

35 ofresiduei 



55 



WO 01/168 



PCT/EP00/08504 



10 



15 



Thus, equation 51 for the weight of a particular rotamer is the partition function that defines 
the probability of having the rotamer z„ from a system of N rotamers. R and T are the gas 
constant and the temperature, respectively- 
Several iterations of field energy calculations (now the weighted average) and 
partitioning are conducted until the set of probabilities is not further modified. It should be 
clear that this/ process is "non-linear, i.e., the quantity to be minimized, MF % depends on 
itself through the adjustable quantities, w. In such circumstances, convergence may be 
achieved through one of several methods of non-linear optimization that will be familiar to 
one skilled in the ait. According to Koehle & Delarue (J. Mol. Biol., 239:249*275, (1994) 
at page 254), convergence is assisted by use of a factor, A. M , set to 0.50. In a preferred 
embodiment of the present invention, in order to break out of dead-locked convergence, X M 
can be set according to the formula (52): 

^-O.SdbO.lRndO (52) 
where Rnd( ) is a random number in the range 0 to 1.0. 



In a preferred embodiment of the present invention, it has been found that smooth 
progress toward the equilibrium state of convergence, is effectively achieved via "simulated 
annealing". This technique is a computer-based method, which simulates the "hcating"of 
20 the protein structure to a high temperature followed by "cooling" it. This is done because 
the high starting temperatures of simulated annealing, e.g., 1000 K, lead to monotonia 
distributions of probabilities of side chain conformations, thus providing a random and 
unbiased starting sample of rotamers. The system is then cooled down to sharpen the 
distributions and optimize the sets of side chain conformations, 



25 



The advantage of the mean field method is that the estimated probabUities are 
similar to frequencies of occurrence, and are coupled to entropy, as shown in the Equation 
(53). 



30 



35 



S ( = -R E^lnw^ (53) 

all rotamers 
t 



wherein S,is the side chain conformational entropy of residue i, R is the gas constant, and 
the set of w represent the probabilities of occurrence of each rotamer. In a preferred 



56 



BNSDOCIO: <WO 01 168U>A2J_> 



WO 01/168 



PCT/EP00/08504 



embodiment, equation (SO) is used to compute entropies of a reference state, as in equation 
(7). 

In a preferred embodiment of the present invention, the mean field approximation is 
used to estimate the weights of all rotamers of the different amino acids in a set of protein 
5 fragments obtained from protein structures in the PDB. In another embodiment, amino 
acids embedded in sample 5-residue extended peptides are employed. From data obtained 
in this way, the reference entropy of the denatured state can be measured. 

The change in entropy of side chains upon protein folding, in both rotamer and sub- 
10 rotamer space can therefore be obtained by a comparison of the entropy of the side chain in 
the native protein, which is then added to the previously determined entropy arising from 
the fixing of the protein backbone (See Section 5.5) to obtain the total change in entropy 
upon folding of the protein. 

15 Mean Field Theory is particularly apposite for the study of amino acid sidechains in 

proteins. In an alternative embodiment, the Mean Field Theory schemes described by Lee 
(J. MoL Biol, (1944) 236:918-939) and Lee and Subbiah, (J. Mol. Biol. (1991) 217:373- 
388), may be employed. Alternative schemes may be employed both for the study of 
proteins and for applications to other macromolecules. In another embodiment, the 

2Q iterative scheme used is Monte Carlo sampling. 



5.11. RE-EVALUATTON OF SOLVATION 
ENERGIES FOR SEQUENCES WITH 

25 

LOW SCORING FUNCTIONS 

In the preferred embodiment of the present invention, if the scoring function of a 
sequence, after being stripped of its solvation term, is below a predetermined energy 
threshold, solvation energies for that sequence are re-evaluated according to two criteria. 
30 Otherwise, the sequence is dropped from consideration. Solvation energies obtained in a 
pairwise manner are removed and re-computed from accessible surface areas derived from 
the optimized configuration of side chains. In the preferred embodiment, the more detailed 
solvation parameters of Eisenberg and McLachlan (1986, Nature 319:199-203), listed in 
Table 2, may be used, though other parameter sets would be adequate. 

35 
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Table 2: Atomic Solvation Parameters 



5 





Atoms 


Radii 
(A) 


Solvation 
Parameters, a, 

(cal/A 1 ) 


Hydrophobic 


C 


1.9 


16 


Atoms 










S 


1.8 


21 


Polar Atoms 


N 


1.7 


-6 




O 


1.4 


-6 


Charged Atoms 


N* 


1.7 


-50 




0 


1.4 


-24 



15 



The accessible surface areas maybe measured using the NSC routine (Eisenhaber et 
al 1995, 7. Comp. Chem. 16:273-284) or another equivalent method. The result is a more 
accurate calculation of the potential energy of the mutant protein. 

20 The solvation energy is assessed according to two properties. As with the previous 

calculation of the energy, the solvation energy of the reference (denatured) state of the 
protein must be considered. This can be calculated for each amino acid by considering the 
solvation energy of a reference state. For this purpose, a reference can be obtained from the 
average solvent-exposure of the amino acid in a 5-residue peptide sequence, as observed in 

25 the Protein Data Bank, but without the context of the surrounding protein structure (except 
for the "capping" residues on the - and C- ends of the sequence. Each residue then behaves 
as if the sequence were a free chain. It is also necessary to consider the environment of the 
residue in the protein- For example, exposed hydrophobic residues should be penalized 
because they may lead to misfolding. Consequently, the solvation energy is calculated by 

30 comparing with residues in all the structures in the PDB. By doing this, it is possible to 
arrive at optimized conformations and sequences for protein-like solvation. 



5.12. OUTPUT OF OPTIMAL SEQUENCE RESULTS 

35 

After the dead-end elimination procedure of the preferred embodiment of the present 
invention, many sequences remain; nevertheless, the subsequent steps (see Mean Field 
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Theory and Refinement, above) ensure that only those conformations and sequences that 
satisfy predetermined energy thresholds finally surface as candidates for the target structure 
The preferred embodiment of the present invention can produce either detailed or limited 
outputs, depending on the size of the sampled sequence space. In one embodiment, the 
5 output is a simple list of sequences and scores (evaluated using the scoring function) that 
can be sorted according to the calculated potential energy so that the lowest energy 
sequences may be readily identified. Li another embodiment, a more complete output 
presents a numerical profile of the energy for each sequence produced. The program is also 
capable of producing a coordinate file (in PDB format) with the structure of the protein 
I0 having a mutated sequence. If mean field sampling is performed, both the PDB-file and 
detailed energy outputs correspond to the combination of most probable rotamers In 
another embodiment, the detailed energy output includes the effective solution score taking 
mto account all candidate rotamers with the weights they were given, and a second detailed 
description of the separate pairwise energy terms resulting of the combination of all possible 
15 sidechamrotamers. If DEE is used for the conformational sampling (without subsequent 
application of MFT afterwards), then the effective solution score corresponds to the GMEC 
where one and only one rotamer is retained for each amino acid side chain; the detailed 
energy file offers nothing else than the separate energy terms which produce the GMEC 
total energy and the PDB coordinate file represents the GMEC model. 



S - 13 - PP NERAIJ7.ATKW TO Ttfn^.pEPTmPS 

Whilst the foregoing has focused specifically upon the types of macromolecules 
25 known as proteins and their building blocks, amino acids, the methods can readily be 
applied by one skilled in the art to other categories of macromoleoules including, but not 
hmted to, those which can be viewed as comprising a fixed structure attached to which are 
freely rotating groups. The alternative embodiments which would be required concern the 
nature of the rotamer library, the choice of reference state, and the property of interest to 
30 optimize. Even though amino acid residues have a simplifying feature in that they consist 
of stde chain and a fixed backbone which enables their conformations to be simply 
expressed as rotamer libraries, other building blocks may also be conveniently modeled by 
one stalled m the art Sugar molecules and nucleotide bases have freely rotating groups 
attached to ring systems, simplifying structural features which would permit the 
35 straightforward construction of conform er libraries. Conformed can be obtained from 
known structures of carbohydrates and nucleic acids respectively or modeled 
computationally. Published molecular mechanics parameters are based on atom-type only 
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and therefore in many cases can be utilized in classes of molecules other than those for 
which they were parameterized. 

The idea of a reference state, although usefully expressed as the denatured form 
when modeling proteins, can be defined differently for other molecules. By analogy with 
5 the alanine pentapeptide, a reference saccharide molecule or small sequence of nucleotide 
bases could be established as a reference structure for carbohydrate or nucleic acid 
modeling, respectively, in a manner similar to the procedure already described. In other 
applications it may be useful to utilize an unsolvated molecule as the reference. 

10 Solubility itself may be a property that can be the subject of investigation and 

optimization with Perla. 

In applications where the property of interest is an interaction between the target 
molecule and some binding molecule, the reference state can be considered to be the 
15 unbound target molecule or a sum of contributions from the unbound target molecule and 
unbound binding molecule. This type of application is likely to be widespread, for example: 
the interaction between DNA and a protein (e.g,, a transcription factor); RNA and a protein; 
the interaction between peptides in solution; the interaction of a polar macromolecule and a 
lipophilic membrane. 

20 



6. EXAMPLE 

Example 1 

Parts of the SH3 domain of a-spectrin, a small globular protein domain with a p-sheet 
architecture, were re-designed with Perla . Nine residues in the buried core of the domain 
were replaced by different hydrophobic amino acids. After the evaluation of template - side 
chain and side chain - side chain interaction energies, the large sequence space was reduced 
by eliminating dead-ending amino acids. For all remaining sequences, the mean field 
approximation was used to set the weights of all possible side chain retainers, providing at 
the same time an estimation of the change (upon folding) in the entropy of the side chains. 
The most probable conformation then served for the computation of the change (upon 
folding) of the protein solvation. A comparison of the best candidates for the replacement of 
the nine wild-type residues indicated that most of them indeed arc constrained by the 
template and do not tolerate mutations, while others undergo conservative changes. Four 
residues that form a surface-exposed turn were similarly re-designed, allowing both 
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nonpolar and polar amino acids. Perla produced a unique sequence, also related to the 
original wild-type sequence. 

Choice of template: The three-dimensional architecture (template), residue numbering and 
wild-type sequence used in the design correspond to the structure presented by Musaccio et 
at, 1992, (pdb accession code lshg; Musacchio, A., Noble. M., Pauptit, R, Wierenga, R 
and Saraste, M. Crystal structure of a Srchomology 3 (SH3) domain. Nature. 359,851-855). 

Operating Parameters: 

1 o The scor »ng function: For the various energy terms, recommended weights were set at 1 .0. 
The force field employed was ECEPP. 

When calculating molecular mechanics energy terms, an interaction cutoff was employed. 
(It is well-established that pairwise interactions between atoms separated by greater than a 
certain distance contribute negligibly.) Here 20 A is found to be large enough to avoid any 
important truncation of the electrostatic interaction energy; van der Waals interactions are 
only taken into account for atom-atom distances smaller than 8 A. 



15 



If the coulomb equation is multiplied by a decaying exponential, the energy reaches zero 
faster than with the alternative scheme and a smaller cut off might be used. Perla selects 

20 the screening scheme according to the value of the screening factor. For large values (i.e. 
more than 1000), a distance-dependent dielectric constant is used. On the contrary, the 
dielectric constant is fixed and the coulomb equation is multiplied by a decaying 
experiential, whose decay factor (1/rJ is set equal to the inverse of the screening factor. 
The dielectric constants, in our example, are 16 and 4, for the solvent-exposed and buried 

25 residues, respectively. Dielectric constants should not be less than half nor more than twice 
those given in this example, otherwise the importance of the electrostatic term would not be 
proper (we had a series of good results using constants of 8 and 4 to measure interaction 
energies between side chains in o-helices). 

30 Modeling sets. In the following examples. Perla takes different amino acid side chains from 
its rotamer library and assembles them op top of nine buried, or four exposed, positions of 
the SH3 domain. The first modeling set (called CORE) consists of v9, Al 1. V23, M25, L31, 
L33. V44. V53 and V58. The second set (called SURFACE) comprises V46. N47, D48 and 
R49. Since all other side chains are kept in the conformation deposited at the Brookhaven 

35 data bank, they constitute the protein template, along with the main chain f the whole 
protein. 
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Amino acids considered. For the CORE modeling set, only nonpolar residues (AVILFW) 
were considered to speed up the sequence sampling, through reducing the total number of 
sequences. Polar and charged amino acids could be avoided since all residues were to be 
fully buried. Since there were 9 residues to design, the total number of sequences was 10 7 . 

? For the SURFACE modeling set, both polar and nonpolar amino acids were considered 
(AVELGDNSTEKRYW); total number of sequences 10 4 * 7 . 

For the CORE set, the amino acid considered have 1, 3, 9, 10, 9 and 12 rotamere, 
respectively, which means that 44 side chains are modeled onto the nine positions, resulting 
10 in 396 constructions. A similar calculation indicates that, with a total of 1400 chain 
confonners, the second design set shows more conformational complexity. 

Solvation of the WT protein, used to fix the threshold for good solvation. The solvation 
energy of the protein target, measured with the sequence and side chain conformations of 
the wild-type protein, is 4,26 kcal mcl" 1 (using the parameters of Eisenberg and McLachlan, 
1986). This should be compared to the solvation energy state, or some other reference, for a 
peptide chain of the same sequence composition. Obviously, the amino acids are more 
exposed to the solvent in the denatured state. Taking as a reference 5-peptides extracted 
from the protein data bank, to mimic the denatured conformations, the solvation energy of 

2Q the wild-type sequence is much higher. The difference in solvation energy between the two 
essential states of the protein folding reaction is thus in favor of the folded three- 
dimensional structure. This should always be expected, since the solvation energy term is 
thought to represent the hydrophobic effect that leads to the compaction of protein chains. 
Besides, taking the average solvation energy calculated in the protein data bank itself, the 

25 difference in solvation energy appears as a penalty. That should be the trend for small 
proteins, for which the buried surface areas are relatively smaller than those of the proteins 
distributed in the data bank. 



Table 3. Original solvation energy (kcal mol' 1 ) (WT sequence with respect to the two possible reference 
states of the solvation potential 



30 



Target 


AM 










Reference 


Difference 




S-pcptide 


32.37 


•28.11 




Protein 


-1.97 


6.23 



35 

In addition, for proteins that have a high proportion of turns, the burial of charged amino 
acids to satisfy some of the main chain hydrogen bond donors is quite common. That is the 
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case of Glu22, which forms a hydrogen bond to the side chain of the amide group of Serl9. 
The impact on the solvation energy term is negative. 

Reconstruction of side chains, estimation of the template-side chain energy term (van 
der Waals). Evidently, the hydrogen bonding term plays no role in the design of the CORE 
set since none of the amino acids considered displays any hydrogen bonding aptitudes, 
while the electrostatics term has a minor contribution that ensues from the uncharged nature 
of the hydrophobic amino acids. The major term for this modeling set is the van der Waals 
term, which Is worth several kcal mol" 1 . The different behaviors for valines at positions 1 1 
and 44 give an interesting hint about the correlation between the molecular mechanics 
potential and the rotaraer preferences found in the protein data bank. Valine, a p-branched 
amino acid, populates preferably the extended conformation, matted by values of \|f above 
100°, where the optimal side chain rotamer is the trans conformer. That matches the case of 
position 1 1 ( <J>--80 and i|r- 120). However, for a value of tj; above 150-160° that coincides 
to the situation at position 44 (<|>=-140 and t|i = 180°, the most abundant rotamer corresponds 
to the gauche - conformer. Clearly, the van der Waals energy term measures these different 
propensities. 



20 Tabic 4. Interaction with the template (kcal mol* 1 ), modeling set CORE 



Ala 9 



He 9 

25 



30 

W 11 gauche - 

trans 

35 

gauche + 



Van der Waals 


Electrostatics 


H-bortding 


-5.78 


0.01 




-5.77 


0.0] 




>50 


-0.03 




>50 


-0,03 




5.27 






-1.24 


0.01 




-13.25 






-13.17 






45.92 


-0.10 




44.76 


-0.11 




-IMS 


-0.05 




-11,21 


-0.05 




42J2 


-0.11 




39.78 


-0.12 
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5 



10 



Val44 


gauche - 


-10.21 


-O.U> 






-10.17 


•0.05 




trans 


0.07 


-0.05 






-0.56 


-0.06 




gauche + 


*7 01 


A AT 

-U.vO 






-3.46 


-0.06 


Trp44 




161.27 


0.03 






98.27 


0.02 






-16.02 


-0.04 






-18.41 


-0.04 



For each angle for each residue the first and second lines represent the energies prior 
and after energy optimisation (achieved through the sampling of subrotamers 
configurations), respectively. Interaction energies of only a few retainers are shown. Zeros 
were replaced with - for clarity. 



Table 5. Interaction with the template (kcal mar 1 ), modeling set SURFACE 

Van der Waats Electrostatics H-bonding 

8.0' 



25 



30 



Lys48 


•3.71 


-0.64 






-355 


-0.77 






-4.29 


•0.47 




Asp 48 


-2.68 


0.67 






-3.08 


0.43 






-3.70 


0.91 




Arg49 


-8.92 


-3.16 


-23S 




-9.77 


-1.06 






-833 


-1.27 


-0.58 


Lys4S 


-3.71 


-0.33 






-355 


•O.40 






429 


-0.24 




Asp 48 


•2.68 


034 
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-3.07 0.22 
-3.71 0.47 



Arg49 -8.96 -1.63 -1.91 

-9.76 -0.54 
-8.33 -0.65 -0.51 



Calculations were run using a dielectric constant of 8 or 16. Only the optimised interaction energies 
arc shown, for a few rotamers. Zeros were replaeed with - for clarity. 



Importance of the reference state. Inspecting the van der Waals interaction energy due to 
1 5 different amino acid side chains, a clear correlation with the size of the chain is found: while 
the small alanine (one heavy atom) scores about six kcal mol* 1 , isolcucine (four heavy 
atoms) contributes approximately 12 and tryptophan (imidazole ring) 18. In the many 
configurations of the denatured state, contacts between the same side chains and the rest of 
the protein certainly exist, and a similar scaling should apply. Hence, since we are interested 
20 only in the energy difference between denatured and folded states, it is necessary to estimate 
and remove amino acid type-dependent reference values from the template-related energy 
term. These are about 3, 5 and 10 kcal mol" 1 for alanine, isoleucine and tryptophan, 
respectively. 

25 Reconstruction of side chains, estimation of the template-side chain energy term (H- 
bonds and electrostatics). To focus on the electrostatic and H-bonding energies, a couple 
of examples from the second modeling set are given (see table 5). Calculations were 
performed with different distance-dependent dielectric constants to review the importance f 
that parameter. Basic amino adds (lysine and arginine) fit better than acidic residues at 

3Q positions 48 and 49. This is a consequence of the overall negative charge of the protein 
template in the proximity of the engineered turn. The amplitude of the electrostatic term for 
residues that are feeing the aqueous solvent is the main point of interest, taking into account 
that water molecules normally screen atomic charges. The electrostatic and hydrogen- 
bonding interaction energies that correspond to the formation of a salt-bridge between some 

35 rotamers of arginine 49 and glutamate 17 of the template, seem to be overestimated when a 
dielectric constant of eight is used. Besides, replacing aspartate 48 by a lysine would 
change the electrostatic term by up t 1 .7 kcal mor\ a value that is very likely to be 
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attenuated by the shielding solvent. Values that are more adequate are obtained with a 
dielectric constant of 16 or even 32. 

Sampling of sub-rotamers to optimize the energy of interaction. In the first example, 
Table 4, the optimization was achieved by averaging the energy of the rotamer found in the* 

5 library and other subrotamers, which combine 5° rotations around the bonds that define the 
X\ and X2 dihedral angles. In some cases, the energy does increase (e.g. the gauche - 
conformer of valine at position 44) due to the overall less favorable interaction of the 
additional subrotamers. For other confonners. slight to significant improvements arc 
obtained (e.g. -0.53 kcal mol" 1 for the gauche + of the same residue and -6.51 kcal mol" 1 

^ for an isoleucine rotamer at position 9). Importantly, large and positive values of the van dex 
Waals term reveal side chain rotamers that are not compatible with the protein template: the 
Gauche confonners of valine at position 1 1 and some conformers of isoleucine, position 9 
or tryptophan, position 44. The repulsions are not removed with 5° rotations of the side 
chains. Some caution must be taken, since side chains can deviate by more than five degrees 

^ from their ideal conformation, as manifested during the analysis of the protein databank 
that led to the generation of the rotamer library. 

The optimization should be given more flexibility; a larger range of subrotamers should be 
sampled. In the following table, some cases are analyzed using two or three 5° steps. The 

20 strong repulsions are not removed, even if the algorithm rotates the valine side chain by 15°. 
The most striking data is that the energy of interaction of the trans and gauche + rotamers, 
with the template, can be significantly improved if alternative conformations we sampled 
within the range of dihedral angles observed in the protein data bank. It is essential to decide 
which scheme of optimization (how many steps, what step size) should be used for the 

25 design of a sequence, but it is extremely difficult to assess the quality and validity of the 
results. Energy changes are important and already in the range of protein stabilities (few 
kcal mol* 5 ) and shall have a serious impact on the final score of each sequence. The third 
optimization scheme should be used with care, Indeed, the description of the sequence 
score by means of template-related energies and side chain pairwise interactions, optimized 

30 separately, might lead to a combination of side chains that fit well when taken by pairs, but 
cannot stand together in the template. The risk for such an occurrence is related to the 
optimization procedure: the larger the rotations, the higher the risk. 

Tabic 6. Optimisation of the interaction with the template (in local kcal mol 1 ) 

35 2x5* ~ _ 

Van der Waals Electrostatics H^anding 
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IS 



Vol 11 


gauche- 


43.28 


-0.13 






trans 


-11.25 


-0.05 


- 




gauche + 


36.90 


-0.13 


- 


Val44 


gauche- 


-10.14 


-0.05 


- 




rrans 


-1.S2 


-0.06 


- 




gauche + 


-3.85 


-0.06 
3x5° 


- 






Van der Waals 


Electrostatics 


H-bonding 


vain 


gauche - 


42.07 


-0.14 






trans 


-11.20 


-0.05 






gauche + 


34.42 


-0.13 




Val44 


gauche - 


-10.10 


-0.05 






trans 


•242 


-0.06 






gauche + 


-4.05 


•0,05 





Main chain entropy loss. When side chain are constructed upon any position, the identity 
of the engineered amino acid is associated to the <f>- jangles to determine the cost of fixing 

20 the backbone, which participates to the entropy change upon protein folding. Although that 
part of the potential acts similarly to secondary structure propensities, only positive energies 
are obtained due to the employed mathematical formulation: good secondary structure 
formers show lower costs (e.g. -0.7 kcal mol" 1 for alanine in a-helical conformation and - 
1 .4 kcal mol" 1 for valine in a P-strand) than breakers (e.g ~1 .4 kcal mol" 1 for glycine in a- 

25 helical conformation and -2.9 kcal mol* 1 in a p-strand). Rarely observed conformations are 
thought to bring tensions into the protein structure, and are marked by high entropy cost, 
whatever the amino acid is (see position 47, following table, and by repulsive van der Waals 
interactions (both terms probably overlap). 

30 Table 7. Interaction with the template (kcal mof* ): entr opy costs an d reference state 

Energy Entropy Cost Reference Sum 

•10.63 -3.07 
-11.92 2.55 5.01 -436 

-10.31 -2.75 

67 



Val46 

35 

Gfy47 



BNSDOCtD: <WO 01 16810A2_I_> 



WO 01/1681 



PCT/EP00/08504 



-0.58 

As* .47 

14.66 
13.89 
49.16 
22.87 
14.26 



4.76 



4.71 



0.84 



5.80 



5.02 

25.16 
24.40 
59.67. 
3337 
24.77 



15 



20 



25 



Elimination of rotamers not compatible with the template. Obviously, confonners with 
strong clashes can be removed from the modeling set. We can deduce from the examples 
shown above that an energy threshold of 0.0 kcal mol" 1 would probably work fine. To 
illustrate the limitation of such an absolute threshold, an example will be drawn from the 
second modeling set, focusing on position 46 and 47. Having a first look at the possible 
conformations of the wild-type valine at position 46, we see that the template-related 
energies, including van der Waals, electrostatics, H-bonding and the cost for fixing the main 
chain dihedral angles, represent favorable interaction energies, even after the reference state 
contribution is removed. There would be no reason at this stage to eliminate one or the other 
rotamer, say that which interacts less favorably with the template (third rotamer, actually the 
gauche + confonner). That same rotamer might be the one that fits better to the target 
protein structure, If it provides an extra few stabilizing kcal mol" 1 through the interaction 
with all other modeled side chains. Consequently, an absolute energy threshold about zero 
would de excellently since all three rotamers would be maintained. 



As shown in the previous table, the template-related energies of an Asn or a Gly residue 
placed at position 47 follow a distinct trend. For glycine, which has no side chain, the 
rotamer built by Perla corresponds to the second alpha proton. It scores few positive kcal 

3° mol* 1 because the molecular mechanics terms (van der Waals, electrostatics and H-bonding) 
are largely overwhelmed by the cost of fixing the main chain dihedral angles. At this 
position, the entropy cost is quite high regardless of the amino acid type (above 4 kcal 
mol* 1 ), due to the 4>-\|f dihedral angles that occupy a less-favored region of the 
Ramachandran plot ((f)* 3 40 and tJj=-100). All energies related to asparagine retainers are 

35 largely positive, mostly because of repulsions with the main chain atoms. If the energy 
threshold for elimination were set up at zero, as previously suggested, neither glycine nor 
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asparagine would be allowed. None of the 20 amino acids would be accepted; Perla could 
not go on with the design. The same absolute threshold simply cannot be used at all 
positions of the target protein. Neither would it be practical to place position-dependent 
thresholds, for several reasons. First, the problem just illustrated would not be completely 
eliminated. Second, there is no simple method to decide what threshold should correspond 
to which position. Third, different thresholds could introduces position-dependent biases, 
for example by allowing more or less amino acids types and conformations, 



Table 8. Elimination of side chains that are the least compatible with the template 



15 



25 





Energy (kcal moT 1 ) 


Weights 


Decision 


Ash 47 










25.16 


0.14 






24.40 


0.50 






59.67 


<0.001 


dismissed 




33.37 


<0.001 


dismissed 




24.77 


0.26 






209.66 


<0.001 


dismissed 




514.07 


<0.001 


dismissed 




403.94 


O.001 


dismissed 




163.95 


<0.001 


dismissed 




223.43 


O.001 


dismissed 



Turning the rotamer-template interaction energies into weights in such a way that 
information is obtained about an amino acid rotamer with respect to the existing alternative 
rotamers (through a partition function), elimination is done without taking care about the 

30 position-dependent energy ranges. Perla has no other choice than keeping the single rotamer 
of glycine (whose weight is one), while the worst performing conformers of asparagine are 
dismissed. The absolute threshold is used, still, as a first barrier (50 kcal mol' 1 ) and a lower 
limit is set for the relative weights (usually 0.001). Both modeling examples were executed 
with the recommended threshold values, sampling subrotamers either with one, two or three 

35 5° steps around the %\ and % 2 dihedral angles of the library elements. Elimination of side 
chain conformers than are not (or least) compatible with the template led t a great 
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reduction of the combinatorial space, therefore, accelerating the stage of side chain - side 
chain pairwise energy calculations. 

Solvation. The presentation of separate intrinsic and pairwise contributions to the solvation 
energy is not of interest- Solvation energies were computed with values of X for core and 
non-core residues, as described above. 

Side chain - side chain interaction energies. These are measured and optimized with 
subrotamers similarly to the way in which the template side chain interaction energies were 
treated. 

Dead-end elimination of sequences. With all possible energy terms stored in memory, 
Perla's task is to set a score for every sequence by summing over all corresponding terms. 
Due to the large number of sequences, and the large number of side chain conformations 
related to each of them, the computational time would be far to long if there was no means 
to skip uninteresting sequences, The dead-end elimination, a mathematical theorem based 
on the pairwise formulation of the scoring function, enables the discovery of side chain 
conformers that cannot belong to the global minimum energy conformation. Such 
conformers can then be ignored. 

Table 9« Dead-end elimination, modeling set CORE 

Number of retainers Amino Acids 



25 



30 



35 



Residue 


Before 


After 


Before 


After 


Position 










9 


9 


5 


(AVILFW) 


AVIL 


11 


6 


3 


U 


AVI 


23 


8 


4 


it 


VIL 


2S 


9 


6 


«« 


AVIL 


31 


15 


5 


*l 


LI 


33 


13 


3 




LI 


44 


8 


3 


At 


AVI 


53 


9 


6 


tt 


AVUJF 


58 


15 


6 




AVLF 




Number of conformations 


Number of sequences 




10" 


10" 


10™ 


10 4J 
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Table 10. Dead-end elimination, modeling set SURFACE 



Number of rotamers Amino Acids 

Residue Before After Before After 

Position 

46 61 1 AVILGDNSTEQKRY V 
10 W 

47 96 1 G 

48 153 1 " S 

49 136 1 K 
j j Number of conformations Number of sequences 

10 w 1 10 17 1 



Mean field approximation, conformational sampling. The mean field approximation sets 
weights for all existing side chain confonners, depending on the pressure maintained by th 
surrounding environment (field). Since that exact environment is itself variable, the 
methodology again costs in an iterative process. Initial weights are established according to 
the interaction of rotamers with all possible environments given at first equal opportunity. 
After computing the weights of the rotamers at every position, the interaction between any 
rotamer with the rest of the modeled side chains is averaged following the weights of their 
own rotamers. With the new field energies, new weights are obtained. The calculation is 
repeated until equilibrium, which is indicated by an insignificant variance of the weights. 
That scheme of cycles if first carried out at high temperature. The system is cooled down 
after each convergence point, until the temperature used to determine the probability 
distribution of the reference state is reached. At that stage the energy score of the sequence 
correspond to a weighted average integrated over all possible side chain conformations. 
Tables below give some details about the score evolution and entropy change during the 
simulated annealing run. The entropy of the side chains, derived from the distribution of 
probabilities, decreases parallel with the temperature decrease. Simultaneously, the fitness 
score improves due to a larger contribution from low energy conformations (e.g, the 
GMEC). The number of cycles necessary to reach a stable set of weights increases, possibly 
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because the energy landscape determined along the multi-dimensional conformation space 
becomes rougher and rougher. 



5 Table 11. Mean field approximation, modeling set CORE (sequence IVKLLVTV with 2520 
conformations) - — — — — — — 



10 



15 



Temperature (K) 


Number of 
cycles 


Absolute entropy 
(kcal mol"') 


Weighted Score 
(kcal mor 1 ) 


1073 


15 


9.04 


-86.76 


973 


23 


7.98 


-87.04 


573 


58 


3.88 


-88.49 . 


473 


67 


2.98 


-88.91 


373 


77 


229 


-89.27 


303 


86 


2.01 


-89.42 



GMEC 
(kcal mol' 1 ) 



-70.77 



20 



25 



1BUK 1*. **w» 

Temperature (K) 


Number of 
cycles 


Absolute entropy 
(kcal mol') 


Weighted Score 
(kcal mol ') 


GMEC 
(kcal mot*) 


1073 


16 


12.60 


5.11 




973 


23 


11.26 


4.94 




873 


30 


9.93 


4.75 




473 


63 


4.56 


3.66 


2.06 


373 


72 


3.26 


3.28 




303 


81 


2.38 


2.99 





30 Thresholds for elimination of sequences. Reconstructing the nine wild-type residues 
corresponding to the CORE modeling set (VAVMLLWV), the sequence-to-structure score 
excluding the solvation term is -48.47 kcal mor 1 . Hence a threshold of -40.0 kcal mol" 1 
would offer a good first selection leveL The solvation term then contributes about the same 
amount of energy as the original X-ray structure, pushing the fitness score below -70 kcal 

35 mol '. F r me second modeling set, the first energy score is positive, due to the strain related 
to the main chain configuration of the turn. The solvation term was slightly improved thanks 
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to the side chain placement carried out by Perla. To design the turn sequences, the two 
energy thresholds can be set to 30.0 and 0,0 kcal mol' 1 , respectively. 



Table 13. Sequence-to-structure fitness score (kcal mol" 1 ) 



Input Target 



10 



Energy 
Solvation term 
Total 



-28,11 



Wild-type Sequence 

Modeling set 

Modeling set CORE SURFACE 
-48.47 24.30 



-28.23 
-76.70 



-3051 
-6.61 



Score of the sequences. Many possible sequences had to be sampled for the CORE design 
j 5 example. When such a case is expected, it is convenient to simplify the output data written 
by Perla, and ask the program to produce a single output file with a list of sequences and 
fitness scores, rather than many detailed energy files. Also, there is no need to have a PDB 
file for every sequence, since these could not be examined. The program can be executed 
afterwards with a reduced set of amino acids or more stringent acceptance threshold levels, 
2Q in order to obtain the fully detailed output. In figure 13, the fitness score is plotted versus 
the ordered list of sequences designed for the CORE 

Designed sequences. Sequences engineered by Perla resemble the wild-type (WT) 
sequence, for both design examples. Out of nine residues, four amino acids were preserved 
in the CORE set (best sequence IVIILLVIV), Only one out of four is maintained for the 
SURFACE example, but the Asp48 to serine and Arg48 to lysine mutations are conservative 
(unique designed sequence VGSK). 

Output Sequences. Sample sequences along with their solution scores as output from the 
program are as follows: 

30 

Table 14 

CORE SURFACE 

VAVMLLWV -76.6 (Wild Type) VNDR -6.6 (Wild Type) 

LVTVLLVIV -81.8 VGSK -28.0 

35 WHLLVIV -8J.8 

IVIILLVW -81.9 
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5 



LHVLLVIV 


-82.0 


IWEIATV 


-82.8 


mVLLVTV 


-82.8 


rvLiLLvrv 


-82.9 


Lvmxvrv 


-83.2 


IVIILLVTV 


-84.4 



10 We note that the solution scores of the best solutions are all somewhat lower than the score 
of the 'Vild type" sequence. Figure 13 shows the distribution of solution scores for 
approximately 1600 considered solutions. 

EXAMPLE 2. Computer-assisted re-design of the SH3 domain and experimental 
characterization of the mutant proteins 

15 



Perla was used to re-design the nonpolar core of the SH3 domain of u-spectrin, (Musacchio, 
A., Noble, M., Pauptit, R., Wierenga, R. and Saraste, M. (1992) Crystal structure of a Sit- 
homology 3 (SID) domain- Nature, 359, 85 1-855), and four solvent-exposed turns. Sequences 
engineered by the design algorithm could be interpreted in terms of packing optimization, high 
secondary structure propensities, and fevorable hydrogen bonding and electrostatics 
interactions. Protein mutants were produced and characterized using circular dichroism, and 
urea-induced equilibrium unfolding was monitored by fluorescence to determine the relative 
protein stabilities. Most mutants do fold to the desired target conformation, and some have 
higher stabilities than the wild-type protein- 



These protein design applications were based on an early version of the computer program that 
30 used only the dead-end elimination principle as a basis for sampling (the technique thus 
resembled that of Dahiyat, B.L and Mayo, S.L. (1996) Protein design automation and Protein 
ScU 5, 89S-903 and Dahiyat, B.L, Gordon, D.B. and Mayo, S.L. (1997)). Automated design 
of the surface positions of protein helices. Protein Sci, 6, 1333-1337. Hence, a unique 
conformation, that of minimal energy, was associated with each sequence, and there was no 
35 side chain entropy term. The optimization of interaction energies was based on conjugate- 
gradient minimisation instead of subro tamers sampling. Solvation was not yet accounted for 
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inthepaixwisedes^^ 

similar to what was described in Example 1 : sequence designs were conducted fixing all side 
chains in the conformation of the X-ray template structure, except for the reaidues to be 
replaced, and similar choices of amino acids were considered. 



Re-design of the buried protein core 

Toreplacethenineresidues buried into the proteincore, we letPerla choose Ala, Val, He, Leu, 
10 Phe or Trp. An example corresponding to this computational problem was presented in 
EXAMPLEl:meo P timalc^ 

withrespecttothewild-typesequenceCVAVMLLVW). Using the early version of Perla, we 
picked up a somewhat different sequence. WLILLVIL, which also contains five mutations. 
Many other sequences were proposed (the dead-end elimination of non-optimal ammo acids 
„ being not totally efficient). As for the example in EXAMPLE 1, amino acid elimination was 
not equal for me nine residue positions (Table 15). Only tryptophan was eliminated at all 
positions,probablybecauseits larger side chain could not fit into the spatially restnctedmternal 
region of this SH3 protein domain. 



20 



t«m» i * Amin o acids participa ting to the most optimal sequences (Core cluster) _ 

Position Wild-type residue Initial choice OpUmal sequences 

9 Val [A,V,I,L,F,W] A.V.LL 



11 Ala 
23 Val 



VAL 
V^,L 
AJ.L 



~ r 25 Met „ rTT 

25 31 Leu " 



33 Leu 



AAL 



44 Val " A.V.LL 

« Val " A.V4p 

58 



Val " Yi^- 



30 



The energy thresholds that help reduce the number of sequences given by the design algorithm 
limitedtheoutputto approximately 3000 sequences (out ofthelO 7 sequences initially possible). 
We decided to construct two protein mutants, relying on the sequence with optimal sequence- 
to-sliucture relationship and on another with less optimal score. This second sequence 
35 (WLLLLAFL) was selected because itprovided one more mutation and contained an aromatic 
residue (Phe 53) whose larger side chain produces a higher impact n the geometrical 
organizationoffcendghbormgsidecha^ 
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aPhe atposition53.au ofthem being correlated with thepresenceof Alain the facing position 
44 When comparing the SH3 domain X-ray structure with the structures modeled by Perla, 
forthese twonew core sequences, we observed that the packing density is increased from wild- 
type to first mutant to second mutant. 



10 



Re-design of the turns 
Fourtumswaedesignedind^ 

acids taken in the following list: A. G, V. I, L, D, E, K, R, H, S, T. N, Q, W and Y (total 
number of sequences: 65536). An example corresponding to the modeling of residues 46 to 49 
(distal loop) was presented in EXAMPLE 1. Using a version of Perla, we obtained a single 
optimal sequence for three of the turns (diverging rum, n-Src and distal loop of the SH3 
domain), while 1 8 candidate sequences were possible for the fourth one (RT-loop). Table 16 
j j presents the designed sequences. 



20 



Table 16. Designed mutants for the four rums. 


Cluster (residues) 


Wild-type sequence 


Designed sequence 


RT-loop (19-22) 




R.S.D.E 


diverging turn (26-29) 


K,K,G,D 


R.H.G.D 


n-Src loop (38-41) 




D,K,D,R 


distal loop (46-49) 


V,N,D,R 


I,G,T,K 



25 Experimental characterization of the protein mutants 

Re-designed proteins were producedby means of recombinant DNA technology and molecular 
biology techniques. Protein expression was not especially high, yet protein yields after 
purification were sufficient to analyze all protein mutants except that corresponding to the 
30 diverging turn (residues 26-29) and n-Src loop (residues 38-41). To check that the proteins 
were correctly folded, far UV CD spectra were recorded. Equilibrium unfolding transitions 
were examined to evaluate the protein stabilities. 



35 Structural characterization via far-UV CD spectroscopy 
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Core mutants. The shape of the far UV CD spectrum of the spectrin SH3 domain is unusual: 
it does not show the typical signal of a A*heet containing protein (a minimum around 217nm) 
hut displays two maxima at approximately 220nm and 235nm, separated by a minimum at 
227nm (Figure 14, A). These specific peaks, attributed to the two consecutive tryptophan 
amino acids (W41-W42), give an indication that the correct tertiary structure is formed; they 
disappear when the protein unfolds (see the high temperature spectrum in Figure 14. A). 
Hence, we observe that the first core mutant is correctly folded (Figure 14. B). The second 
mutant, however, is nearly unfolded at 298 K but folds properly at lower temperature, attestmg 
that it is less stable than the WT SH3 domain (Figure 14, C). 

10 Turn mutants. Similarly, the three designed mutants we were able to produce for the turn 
clusters (RT-loop, diverging turn and distal loop) are correctly folded (Figure 15). 



1S Urea-induced equilibrium unfolding 
Toconipareprotemstabili^ 

recording the change in emitted fluorescence. The unfolding curves were fitted with Equation 
(54) (based on the linear extrapolation method; Pace, C.N. (1986) Detennination and analysis 
of urea and guanidine hydrochloride denaturation curves Methods Enzymol, 131. 266-280; 
20 Santoro,M.M.andBolen.D.W.(1988)UiuW^ 

extrapolation method. 1. Unfolding of phenylmethanesulfonyl alpha-chymotrypsin using 
different denaturants. Biochemistry, 27, 8063-8068), which corresponds to a two-state 
unfolding mechanism but has been improved to include a special dependence of the 
fluorescenceof SH3 proteins on urea concentration. In equation (54), the parameters of interest 
are the slope in the transition region ( m u _ f ) and the folding free energy in the absence of 
denaturant ( AG^ ). Figures 16 and 17 show the unfolding curves and the parameters 
obtainedfrorodatafittingtoEquation(54) are summarized in Table. Data fitting was improved 
30 defining explicitly the slope in the transition region to -0.93 (the average of the values 
corresponding to the wild-type protein), as observed by the smaller difference between folding 
free energies obtained from separate experiments. 



25 



35 ^obs = 



(F f + m f [D\Y (F u + m v [D))**p\~ (" A <W * [Wg*! ^)/ RT ] 
1 + exp[- (- A C?2/ + rn[D] + cl[D] 2 )/rt] 
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(54) 



Equation (54) is an expression of the observed fluorescence (F ob ) as a functi n of the 
denaturant concentration ([D]). for proteins that have a two-state unfolding transition. The 
linear dependence on denaturant concentration, of the fluorescence of the folded (f? and 

5 unfolded CFJ states is taken into account, with slopes m, and m. respectively. The protein 
stability or folding free in energy in absence of denaturant is &G*° and W M / is the slope 
in the transition region. The quadratic term, particular to the SH3 domain of a-spectrin Prieto, 
j W ilma^M.,Jimcne Z ,M^ 

10 i„' protein folding and stability: introducing a helical tendency in the all beta-sheet alpha- 
spectrin SH3 domain. JMol Biol. 268, 760-778; Viguera, A.R.. Serrano, L. and Wilmanns, M. 
(1996) Different folding transition states may result in the same native structure. Nat Struct 
Biol, 3, 874-880, is added to improve data fitting (a - 0.008925). R and T are the gas constant 
and temperature, respectively. 

15 

Core mutants. The unfolding curves in Figure 1 6 clearly indicate that the first core mutant is 
. m orestablewanmewUd-typeSH3doma^ Data 
fittmgtoEquatton(54)producedparametersmagreementwimpublisheddate 

20 protein (Table 17). For the first core mutant, stability is increased by 0.7 kcal mol" . For the 
second mutant, the fitting gave an indication of the amount of folding free energy lost (at least 
2.5 kcal mor'), but the absence of the initial part of the curve renders the modeled parameters 
quite imprecise. Slopes in the unfolding transition regions are similar, which suggests that the 
hydrophobic surface areas exposed to solvent (in the unfolded and folded states) have not 

25 changed significantly. 



Turn mutants. Then-Src loop proteinmutant was not successfully purified and the diverging 
turn mutant was not purified in sufficient quantities, hence, only two of our designed protems 
forme turn clusters were characterized (Figure 17). One ofthem, the RT-loop mutant, is just 
sUghuymorestablemanmewild-typeSHSd^rnain. On the other hand, modifying the distal 
loop residues resulted in a considerably larger stabilization. Slopes in the unfolding transition 
regions were again similar to that of the wild-type protein domain. 
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Table 17. Folding free energies of the SH3 domain from ar-spectrin and mutants, and 
prediction by Perla of the free energy changes due to the sequence mutations. 



5 



10 



Protein 




'J 


AG*?/ 


AG^ b 

IT-*/ 


AG?°< 


AG*? - 


AAE PlflB 


Distal loop 


-0.94 


-0.81 


-3.86 


-3.54 


-3.82 


-3.85 


-8.95 


Core 1* 


-0.89 


-0.96 


-3.49 


-3.66 


-3.55 


-3.57 


-7.03 


RT-loop 


-0.85 


-0.77 


-2.81 


-2.58 


-3.06 


-3.10 


-0.83 


Wild-type 


-0.88 


-0.98 


-2.72 


-3.02 


-2.84 


-2.85 




Core 2 nd 


ND 


-0.96 


ND 


-0.24 


-0.0 f 


-0.02 


-4.45 


diverging turn 


insufficient protein yield after expression and purification 


1.8 


n-Src loop 


insufficient protein yield after expression and purification 


2.48 



Parameters obtained from data fitting to Equation 54 and °* d fixing the slope in the 
transition region to -0.93 ; * c andb,d first and second experiments. e Difference in scoring energy 
(mutant - WT). The current version of Perla was used to model the wild-type and designed 
sequences, building the side chains of all residues and with 0.5 weights for all molecular 
mechanics energy terms and entropy changes, distance-dependent dielectric constants of 8 
and 4 for solvent-exposed and buried regions of the protein, 2 steps of 5° rotations to generate 
subrotamers and a sampling temperature of 303K, (ND) Not determined (correct data fitting 
not possible). 



25 EXPERIMENTAL PROTOCOLS, METHODS 

Gene construction and cloning. Nucleic acid sequences were obtained by reverse-translating 
(using the preferred codon usage of£. coli) the protein sequences, either wild-type SH3 domain 
of a-spectrin or the designed mutants. Full-length genes were engineered, thus including the 

3Q region coding for the N-tenninal first 5 residues (MDETG, including the initial methionine), 
which are not observed in the X-ray structures of previously characterized mutants and wild* 
type. All DNA sequences were built using a polymerase chain reaction (PCR) method with 
oligonucleotides synthesized by the EMBL DNA service. Two central oligonucleotides were 
annealed together and polymerization was accomplished to obtain a double-stranded DNA 

35 fragment, which was further elongated and amplified using two external primers. Two stop 
codons (TAG) were introduced at the 3' end of the DNA sequence. Ncol (CCATGG) and 
Hindm (AAGCTT) restriction sites were designed at the 5* and 3 1 termini of the produced 
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DNA, respectively, to allow cloning into pBAT-4 (Paranen. J., Rikkonen. M, Hyvbnen, M. 
andlOlariainen, L. (1996) T7 vectors with modified T71ac promoter for expression fprotems 
in Escherichia coli. Analytical Biochemistry. 236, 371-373) after digestion with the 
corresponding restriction enzymes. Electro.compctentorchemically.competentXL-1 Blue* 
« coli ceUs were transformed with the engineered plasmid vectors, and the cells were grown on 
L-brothplatc containing ampicilin. Positive clones were selected after screeningbyPCR (with 
primers complementary to the vector at the 5' and 3' sides of the cloning site) for bactenal 
templates that generate a DNA fragment of the expeeted size (about 200 base pairs). Gene 
sequences were confirmed by chemical sequencing of the vector cloning site, performed 
0 following the dideoxy-mediated chain termination method of Sanger, (Sanger, F., Nicklen, S. 
and Coulson. AJL (1977) DNA sequencing with chain-terminating inhibitors. Proc Natl Acad 
Sci U SA, 74, 5463-5467) after purification of the vectors from these positive clones. 



15 Protein expression, extraction and purification. Eleetro-competent or chemically-competent 
E coli BL21 (DE3) cells were transformed with the purified expression vectors and grown at 
37*C (61 scale) from a single colony in L-broth medium containing 50mg V ampicilin until the 
culture reached an optical density of -0.6 at 600nm. Isopropyl-b-D-thiogalactopyranoside 
(TPTG) was then added to a final concentration of 40mg l" to induce the expression of the 
20 engineered protein. Cells were harvested three hours later by centrifugation (about 3300 RPM 
for half an hour) and resuspended in lOmM sodium citrate P H 3.5 and lOOmMNaCl, lysedby 
sonication and ultracentrifuged at 37500 RPM for 2h. The resulting pellets (containing 
insoluble proteins) were resuspended in 6M urea, lOmM sodium citrate pH 3.5 and 100 mM 
NaCl, while polyemyleneimime (PET) was added to the soluble fractions (containing soluble 
25 proteins) to precipitate DNA fragments. Both fractions were submitted to a second run of 
ultracentrifugation. Electrophoresis on sodium dodecyl sulfate-poiyacrylamide gels (SDS- 
PAGE) was used to determine which fractions contained the overexpressed proteins. While 
proteins from the insoluble fractions were directly purified, proteins present in the soluble 
fractions were first precipitated with ammonium sulfate in two steps. Ammonium sulfate was 
30 added to the solution to reach a concentration of 30% and the sample was ultracentrifuged (at 
25000 RPM for half an hour). Some additional amount of ammonium sulfate was then added 
to thesolublepartto increase the concentration to 70% and a new centrifugation run performed. 
Pellets were as before resuspended in 6M urea, lOmM sodium citrate pH 3.5 and 100 mM 
NaCl, and the suspensions submitted to ultracentrifugation. Proteins purification was carried 
35 outbyexclusionchromatographym^ 

3.5 and 100 mM NaCl) on aHiLoad 26/60 Superdex 75 column. Fractions with pure protein 
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weredetectedwithSDS-PAGE.pooledtogdheranddiluted 10times(in lOmM sodium citrate 
pH 3 S and 100 mM NaCl) to allow refolding, and were afterwards concentrated by 
ultrafiltration. The purification was repeated once and purity checked in polyacrylamide gels 
(purity estimated to be >95%). After refolding and concentration, protem solutions were 
dialyzed against water at P H3.5. Purity and protein identity were checked by mass 
spectroscopy by the EMBL peptide & protein service. About 10-20mg pure proteins were 
obtained. 



10 Determination of protein concentration. Protein concentrations were detemuned usmg the 
me thodofCHu,S.C.andvonffi^^ 

ftom amino acid sequence data. AnalBiochem, 182, 319-326. Absorbance (A 2S0 ) of tyrosme 
andtryptophanresidueswa^ 

The folded samples were buffered aqueous solutions of the protein and the unfolded samples 
15 similar solutions containing additionally 6M guanidinium hydrochloride (GdnHCl). The 
following equations (55 and 56) were used to determine the extinction coefficient at 280nm 
(^andm e polypeptideconcentration,respectively(wim/theceU 



20 



folded .folded •. 

folded ( -i\ _ Jm = - 2 *? M A N Tyr \2W+ N Trp S690\ (55) 

^280 [M cm )- .unfolded E 280 unfolded [ Tyr irp j 

4 2 80 280 



25 



Far-UV circular dichroism. CD spectra were recorded on a Jasco-7 10 instrument calibrated 
using D-10-camphorsulfonic acid. Measurements were made every O.lnm, with a response 
timeoflsandabandwidthdflIun,atascanspeedof5()Dmn 1 in- , . Protein concentrations were 

30 about lOCIinM; sample* werenotb^ 

temperature was set as indicated in the figures. CD spectra were recorded using a cuvette with 
aO^mmpath. Spectra shown in the text are the average of20 scans, which were corrected for 
the baseline signal. 

35 
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Urea-indaced quilibruim denaturation. Urea titrations were realized using a Jasco-710 
instrument equipped with an automated titration system. A 2ml sample (about 15mM protein 
in 50mM sodium citrate pH3.5 and 50mM sodium chloride) was placed in a cuvette with a path 
length of 1cm. Two software-controlled syringes were used to replace in a stepwise manner 

5 a fix volume of the sample with a urea-containing protein solution (15mM protein in 50mM 
sodium citrate pH3.5. 50mM sodium chloride, an at least 8M urea). Thus, the urea 
concentration was increased while the buffer and protein concentrations were kept constant. 
20 urea injections of 50ml were performed first and followed by 15 injections of 100ml. in 
ordertoperformmeasxn-ementsformeaconcentrationsbetweenOMand6M(meprotocolbemg 

10 limitedbythe2.5mlvolumeofthesyringes). Aftereacbureainjection.the sample was allowed 
to equilibrate for 2 minutes, and the total fluorescence above 340nm was recorded (excitation 
wavelength was 280nm). Constant stirring of the solution was used to facilitate sample 
equilibration. Temperature was maintained at 298K. For the urea-containing protein solution 
and the final protein sample, urea concentrations were determined from refractive index 

15 measurements as indicated by Equation 57 (where /Wis the difference between the refractive 
index of the denaturant-containing sample and the refractive index of the denaturant-free buffer, 
Pace, C.N. (1 986) Determination and analysis of urea and guanidine hydrochloride denaturation 
curves. Methods Enzymol, 131. 266-280. Denaturant concentrations for each measurement 
were calculated according to the experiment protocol and the urea concentration of the urea- 

20 containingprotein solution used for urea additions. The difference between the measured final 
denaturant concentration and the concentration expected due to the unfolding protocol was less 
than 0.2M. The stability of pH was also checked and confirmed at the end of the experiment. 

[ tt rea](A0= H7.66(AiV)+29.753(AJ\0 2 + 185i6(A7v0 3 (57) 

25 7. RF.FFRF.NCES CITED 

All references cited herein are incorporated herein by reference in their entirety and 
for all purposes to the same extent as if each individual publication or patent or patent 
application was specifically and individually indicated to be incorporated by reference in its 
30 entirety for all purposes. 

Many modifications and variations of this invention can be made without departing 
from its spirit and scope, as will be apparent to those skilled in the art. The specific 
embodiments described herein are offered by way of example only, and the invention is to 
35 be limited only by the terms of the appended claims, along with the full scope of equivalents 
to which such claims are entitled. 
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What is Claimed is: 

I. Ameth d for choosing a set of substitute building blocks for a set of positions in 
a target macromoiecule according to whether a calculated solution score is lower than a 
threshold value, the method comprising: 

5 

(a) specifying at least one substitute building block for each position in said set of 
positions to produce a specified set of substitute building blocks; 

(b) for each said substitute building block, 

10 i) determining at least one candidate conformer; 

ii) substituting coordinates of each said candidate conformer or portion 
thereof for coordinates of the building block or portion thereof at said 
position in an atomic structure of said target macromoiecule; and 

1 J (c) minimizing the value of a calculated energy term by adjusting the geometry of 

each said candidate conformer or portion thereof in order to obtain a solution structure; 

(d) calculating a solution score for said solution structure, wherein said solution 
score comprises an entropic term; and 

20 

(e) choosing said specified set of substitute building blocks if said calculated 
solution score is lower than a threshold value. 



25 2. The method of Claim 1 wherein said macromoiecule is a peptide or protein; said 

building blocks are amino acid residues; and each candidate conformer is a side chain 
rotamer selected from a plurality of side chain rotamers. 



30 

3. The method of Claim 2 wherein said calculated solution score comprises a 
difference between a first value corresponding to said solution structure and a second value 
corresponding to a reference structure. 



35 
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4. The method of Claim 3, wherein said first value corresponding to said solution 
structure accounts for interactions between said side chain rotamer and said atomic 
structure, and a sum of interactions between all pairs of all possible side chain rotamers. 

5 

5. The method of Claim 3, wherein said reference structure is a denatured state of 
said solution structure. 



6. The method of Claim 4, further comprising a step of rejecting a side chain 
rotamer when the value of said interactions between said side chain rotamer and said atomic 
structure is greater than a threshold value. 



7. The method of Claim 2, wherein the dihedral angles of said side chain rotamers 
are optimized in step (c). 

20 8. The method of Claim 2, wherein the positions of all main chain atoms of said 

atomic structure, and the positions of all atoms in amino acid side chains that are not 
included in said set of substitute building blocks are held fixed in said atomic structure. 



9. The method of Claim 7, wherein the positions of all atoms in amino acid side 
chains on residues that are not at said set of positions are allowed to vary whilst the dihedral 
angles of said rotamer are optimized. 



° 10. The method of Claim 7, wherein the positions of all main chain atoms of said 

atomic structure are allowed to vary whilst the dihedral angles of said rotamer are 
optimized. 



35 

1 1« The method of Claim 2, wherein said plurality of side chain rotamers is a library 
of predetermined rotamer conformations. 
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10 



12. The method of Claim 2, wherein said plurality of side chain rotamers is derived 
from a continuous distribution of conformations. 



13. The method of Claim 1 , wherein 

said atomic structure includes a representation of the building blocks at each position 
in said set of positions; and 

said atomic structure was determined by a method selected from the group 
consisting of x-ray crystallography, nuclear magnetic resonance spectroscopy, electron 
microscopy, homology modeling, and ab initio molecular modeling. 



j S 14, The method of Claim 1, wherein said atomic structure is an X-ray crystal 

structure of a portion of said macromolecule that comprises said building blocks at each 
position. 



20 15. The method of Claim 14, wherein said X-ray crystal structure was determined at 

a resolution of less than 4.0 Angstroms. 



16. The method of Claim 2, wherein said calculated solution score is calculated 
^ using an empirical scoring function. 



17. The method of Claim 16, wherein said empirical scoring function is a sum of 
3Q energy terms, comprising a template energy of said atomic structure held in a fixed 

geometry, an intrinsic energy of interaction between a candidate side chain retainer and said 
atomic structure held in a fixed geometry and a pairwise energy of interaction between 
possible pairs of side chain rotamers in said substitute set of building blocks. 



35 
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18. The method of Claim 17 wherein said intrinsic energy of interaction is 
computed as either £ focd 5tnjcm# best ride ch«n(/> or X W i,r^fixed structure, side ch«n(i,0 

ro tamers 
r 



5 

wherein JJ w toucnire . <We " an energy of interaction between the atomic structure 

held in a fixed geometry and rotamer r of side chain i, and w„ is a weighting factor. 



10 



15 



19. The method of Claim 17, wherein said pairwise energy of interaction is 
computed as either tide ( i)f best tide ehein(7) 01 

J] X W.V^V^d. eheiu(/.r). .Me eh«i B (y.x) Wherein side ^ side >s 

ro tamers r rotaoien f 
of residue I of residue > 

the energy of interaction between side chain rotamer j and side chain rotamer/, and w u and 
Wjj are weights. 



20 



20. The method of Claim 17, wherein the template energy of said atomic structure 
held in a fixed geometry comprises at least one energy term selected from the group 
consisting of a molecular mechanics potential, a solvation energy, an empirical penalty 
function, and an entropic contribution. > 



25 



30 



21. The method of Claim 20, wherein the template energy of said atomic structure 
held in a fixed geometry comprises a sum of terms whose coefficients are individually 
adjustable weighting factors. 



22. The method of Claim 20, wherein said molecular mechanics potential comprises 
at least one energy term selected from the group consisting of bond length vibrations, bond 
angle bends, the hydrogen bond energy between pairs of hydrogen bond donor and acceptor 
atoms, an lectrostatic interaction energy between pairs of charged atoms, and a van der 
Waals interaction energy between pairs of non-bonded atoms in said atomic structure. 

86 



BNSOOCID: <WQ 01 16810A2J_> 



WO01/ 



PCT/EP00/08504 



23. The method of Claim 22, wherein the van der Waals energy term is expressed 



as 



E' I 

vdw ftonbondcd/y 



[A Ml 

12 ~ 6 



wherein the sum runs over all possible non-bonded atom pairs i and y from said atomic 
structure held in a fixed geometry. 



10 



24. The method of Claim 22, wherein the hydrogen bonding energy between pairs 
of hydrogen bond donor and acceptor atoms is expressed as 



15 



r _ y f Aha &ha 

HS — H- bonded H,A\ 12 fif 



20 



wherein the sum runs over all hydrogen bonds in said atomic structure held in a fixed 
geometry and atoms i and/ are donor and acceptor atoms participating in each of said 
hydrogen bonds. 



25 25. The method of Claim 22, wherein said electrostatic interaction energy between 

pairs of charged atoms is expressed as 



£elec= Z, 



30 ChWBCS U 4x£oe r nj 

wherein the sum runs over all pairs of charged atoms i, and; in said atomic structure held in 
a fixed geometry whose respective charges arc q, and q r 

35 
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26. The method of Claim 20, wherein said entropic contribution comprises at least 
one term selected from the group consisting of a main chain entropy term, a side chain 
rotation entropy term and a side chain vibration entropy term, 

27. The method of Claim 26 wherein said main chain entropy term is 
calculated as 



10 >r amino acidr 



Z_ amino 

subspaccs fVznxr 

entropy V 1.. clo5c »* r ' 

-*Wix, chain *V** £j m ^amino acid/ 



all residues i iV tll^ 

wherein w tnt ?*l ■ is a coefficient, T is temperature, and AT is the number of amino acids of 

main chain 

a particular type in a database of Jcnown protein structures that are found within a specific 
range of <j>, uy angles. 



20 28. The method of Claim 26 wherein said side chain rotation entropy term is 

calculated as 



25 ^ W sidTSiin^ryi ^ 



all residues i 



all fo lamer* r 
V of residue i 



target 
structure 



alt retamcrer 
^ of residue / 



reference 
structure/ 



wherein wS !n is a coefficient, T phya is a temperature, and w r is obtained from a partition 



30 wuwwu ^ i<lo chain 

function. 
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29. The method of Claim 26, wherein said side chain vibration entropy term is 



calculated as: 
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-w. 



vibration 



ode chain 



T p>v> Z Z *V 

all residues /all rotamcrsr 
of residue / 



10 



all sub-ro tamers * 
of ro tamer r 



target 
structure 



all sub- retainers * 
of retainer r d 



reference 
structure / 



wherein W*£%* n is a coefficient, w, is a weight, and w, is obtained from 
function. 



a partition 



15 



30. The method of Claim 20, wherein said solvation 



energy is calculated as 



20 



lUOO 



" 8B,/ ^°°»' 7 rrfcrmcc 



J rrfcrmcc 
structure / 



25 

wherein o, is a parameter specific to atom /. is a coefficient, and ASA, is an 

accessible surface area of atom i and atom / is in said atomic structure. 

30 

31. The method of Claim 30. wherein said parameter specific to atom i reflects the 
properties of a solvent selected from the group consisting of water and an organic solvent. 

35 32. Ihe method of Claim 31 wherein the organic solvent is selected from the group 

consistmg of methanol, methylamine, and dimethyl sulphoxide. 
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33. The mcth d of Claim 20, wherein said empirical penalty function is calculated 

as 

^ ill residues / 

wherein add t is a term representing a probability of occurrence of amino acid i in 
nature, T gtat is temperature, and w" a ' is a coefficient 

10 

34. The method of Claim 5, wherein said reference structure comprises said side 
chain retainer substituted for a side chain in an alanine based penta-peptide. 



15 

35. The method of Claim 5 wherein said reference structure comprises said side 
chain retainer embedded in a fragment of protein taken from an atomic structure of a 
naturally occurring protein or an ensemble of fragments of protein, the populations of which 
are determined either from populations in the naturally occurring proteins or from 
20 computations establishing the potential energy of each fragment and integrating them into a 
partition function. 



36. The method of Claim 5, wherein said reference structure is a denatured state of 
2 * said atomic structure and said side chain vibration entropy term in said reference structure is 
modeled as a uniform distribution of sub-rotamer conformations. 



37. The method of Claim 18, wherein said intrinsic energy of interaction comprises 
at least one energy term selected from the group consisting of a molecular mechanics energy 
term, a solvation energy term, and an entropic contribution. 



38. The method of Claim 37, wherein said intrinsic energy of interaction comprises 
a sum of terms whose coefficients are individually adjustable weighting factors. 
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39. The method of Claim 37, wherein said molecular mechanics energy term 
comprises at least one term selected from the group consisting of the van der Waals energy 
between pairs of non-bonded atoms, the hydrogen bond energy between pairs of hydrogen 
bond donor and acceptor atoms, and the electrostatic interaction energy between pairs of 
^ charged atoms. 



40. The method of Claim 39, wherein the van der Waals energy between pairs of 
non-bonded atoms is calculated as 



10 



Evdw 5> nooboaded/7 



f An Bo 



6 



wherein i and ; represent all possible atom pains comprising atoms i of said atomic structure 
1 5 held in a fixed geometry and atoms j of said side chain rotamer. 



41. The method of Claim 39, wherein the hydrogen bonding energy is calculated as 



20 



Era- ^H-Wdtd H,A 



f \ 

Aha Bha 

„12 10 

r ij r v J 



25 



wherein the sum runs over all hydrogen bonds between said atomic structure held in a fixed 
geometry and said side chain rotamer, and atoms i mdj are respectively donor and acceptor 
atoms participating in each of said hydrogen bonds. 



30 



42. The method of Claim 39, wherein the electrostatic energy between pairs of 
charged atoms is calculated as 



35 



£elec= ^charges ij Aff £Q£rrij 
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10 



wherein the sum runs over all pairs of charged atoms such that atom i is found in said 
atomic structure and atom j is found in said side chain rotamer and whose respective charges 
are q, and q } . 

43. The method of Claim 37, wherein said eutopic contribution comprises at least 
one term selected from the group consisting of a main chain entropy term and a side chain 
vibration entropy tenn. 

44. The method of Claim 43. wherein said main chain entropy term is calculated as 



IS 



sobspaoes ^K20*j»* 
dose to residue fyr 



residue 

w NT* 



residue 



20 wherein H^Wn is a coefficient, w^is obtained from the partition function, is 
physiological temperature, and Nis a number of amino acids of a particular type found in a 
given range of main chain dihedral angles. 



25 



45. The method of Claim 43, wherein said side chain vibration entropy tenn is 
calculated as 

f \ residue 

.vibration ~~ 



30 



vibration <r» 
~ w side chain *phyt 



V gub-rotamen s J 



r '^reference 



target 
structure 



35 wherein vv^^? is a coefficient, and w, is obtained from a partition functi n. 
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46. The method of Claim 37, wherein said solvation term is calculated as 



^solvation 



of side chain 

wherein tf° ]yvioa is a coefficient, and ASA t is the accessible solvent area of atom L 



^ atoms / of side chain 



47. The method of Claim 19, wherein said paiiwise energy of interaction comprises 
10 at least one energy term selected from the group consisting of a molecular mechanics term, a 
solvation energy term, a penalty term, and an entropic contribution. 



48. The method of Claim 47, wherein said pairwise energy of interaction comprises 
1 * a sum of terms whose coefficients are individually adjustable weighting factors. 



49. The method of Claim 47, wherein said molecular mechanics term comprises at 
2Q least one term selected from the group consisting of a van dcr Waals energy term, a 
hydrogen bond energy term and an electrostatic interaction energy term. 



50. The method of Claim 49, wherein the van der Waals energy term is calculated 

25 as 



' nonbcmfcd/y 

30 



'A, Bj\ 

12 ~ 6 

V fa Tij J 



wherein the summation runs over all possible atom pairs comprising atoms i from the first 
side chain ro tamer of one of said pairs of side chain ro tamers and atoms j from the second 
side chain rotamer of one of said pairs of side chain rotamers. 



35 



93 



BNSDOCID: <WO 0116810A2J_> 



WO 01/168 



PCT/EP00/08504 



51. The method of Claim 49, wherein the hydrogen bonding energy term is 
calculated as 



Ehb~ ^H-bonded HA 



Aha Bha 



V r ij r * J 



wherein the summation runs over all possible hydrogen bonds between atom pairs 
comprising atoms i from the first side chain rotamer of one of said pairs of side chain 
10 rotamers and atoms j from the second side chain rotamer of one of said pairs of side chain 



ro tamers. 



52. The method, of Claim 49, wherein the electrostatic interaction energy term is 
calculated as 



20 



^elec- ^charges ij^ £Q£rnj 



wherein i and j represent all possible pairs of charged atoms comprising atoms i, with 
25 charge q„ from the first side chain rotamer of one of said pairs of side chain rotamers and 
atoms;, with charge q p from the second side chain rotamer of one of said pairs of side chain 
rotamers ■ 



30 53. The method of Claim 47 wherein said entropic contribution comprises a side 

chain vibration entropy term. 

54. The method of Claim 53, wherein 6aid side chain vibrauon-entropy term is 
35 calculated according to an energy-based distribution of sub-rotamers. 
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55. The method f Claim 53 wherein said side chain vibrational entropy tennis 
calculated as: 



vibration *rt * 
' W Pairwise *phys' i 



-R I w An lnw 



sub-rotamets A and B t 
of retainer pair X£ 



rotamerpaiMB 
in target structure 



10 



of no tamer /4 



only rotameM 
in target structure 



15 



-R 5>A,liiWir, 



sub-roumersfl, 
of retainer £ 



only retainer £ 
in target structure/ 



20 



56. The method of Claim 55, wherein the weights, w„ of the sub-rotamere in said 
side chain vibration entropy term are obtained by means of a partition function. 



57. The method of Claim 47, wherein said solvation energy term is calculated as a 
difference between an accessible surface area of the side chain in the reference state and the 
25 measured accessible surface area of the side chain conformer substituted into the atomic 
structure and the reference structure is a denatured form of the macromolecule. 



30 



58. The method of Claim 47, wherein said solvation term is calculated 



aorvatioa y J i^^SJ^' 

taultf rdtfA <r » V 1 ^ 'l« 1*50 ««uc*m / 
«f kMn fair A3 



as 



35 
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wherein (ASA,) «Jy residue a«b is a measured accessible surface area of atom i in each side 

V nj in taip* structure 

chain confonner substituted, separately, into the atomic structure, and {ASA^^mg^ 

5 is the measured accessible surface area of atom i in each side chain substituted together into 
the target atomic structure. 



10 



15 



25 



35 



59. The method of Claim 47, wherein said penalty term is calculated as 



Mat 

amino *cid i 

all residues/ 



wherein is a coefficient, and P is a probability of occurrence of a pair of amino acid 

residues whose side chain rotamers are contributing to the pairwise interaction energy. 



20 



60. The method of Claim 17, wherein said intrinsic energy of interaction comprises 
a molecular mechanics term and said pairwise energy of interaction comprises a molecular 
mechanics term, and wherein the dihedral angles of said side chain rotamer are optimized in 
step (c) by minimizing both the molecular mechanics term of said intrinsic energy of 
interaction, and said molecular mechanics term of said pairwise energy of interaction. 



61. The method of Claim 60, wherein the dihedral angles of said side chain rotamer 
30 are optimized in the space of internal rotations of said rotamers by an algorithm selected 

from the group consisting of least squares, steepest descents, quasi-Newtonian, and 
simulated annealing. 

62. The method of Claim 60, wherein the dihedral angles of said side chain rotamer 
are optimized by averaging over the values measured by sampling sub-retainer 
configurations of the side chain rotamers, generating said sub-rotamer configurations by 
sampling a range of dihedral angles by stepping each dihedral angle in said side chain 
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rotamcr by a predetermined step size for a number of steps and, selecting said sub-rotamers 
whose contribution to the calculated solution score is highest. 



5 63. The method of Claim 62, wherein the predetermined step size and the number of 

steps sampled is determined by an amino acid type of the side chain rotamer. 

64. The method of Claim 6, wherein the side chain rotamers which are not rejected 
^ and that have a probability smaller than a predetermined probability threshold are rejected. 



65. The method of Claim 64, wherein said probability is calculated from a partition 
function over all possible rotamers of a particular residue type, using their energy of 
interaction with a fixed portion of the atomic structure wherein said energy of interaction 
comprises at least one energy term selected from the group consisting of a molecular 
mechanics energy term, a solvation energy term, and an entropic contribution. 



20 

66, The method of Claim 2, wherein each candidate confonner specified in step (b) 
is a D-enantiorner selected from the group consisting of glycine, alanine, valine, leucine, 
isoleucine, glutamic acid, aspartic acid, asparagine, glutamine, proline, phenylalanine, 
tyrosine, serine, threonine, lysine, arginine, histidine, cysteine, tryptophan, and methionine. 

25 



67. The method of Claim 2, wherein each candidate confonner specified in step (b) 
is an L-enantiomer selected from the group consisting of: glycine, alanine, valine, leucine, 
isoleucine, glutamic acid, aspartic acid, asparagine, glutamine, proline, phenylalanine, 
tyrosine, serine, threonine, lysine, arginine, histidine, cysteine, tryptophan, and methionine. 



68. The method of Claim 2, wherein each candidate confonner specified in step (b) 
is selected from the group consisting of the L~ and D-enantiomers of amino acids including, 
but not limited to, norvaline, beta-alanine, and tartrine. 
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69, The method f Claim 11, wherein said library of predetermined rotamer 
conformations is constructed by a method comprising: 

tabulating, for each of the twenty naturally occurring amino acids, a 
statistical distribution of observed amino acid side chain dihedral angles in a set of 
ciystallographically determined protein structures; 

determining a Gaussian distribution of observed amino acid side chain 
dihedral angles tabulated in said tabulating step; and 

constructing amino acid side chain rotamers for each of the twenty naturally 
occurring amino acids using all combinations of Gaussian peaks around each amino 
acid side chain dihedral angle tabulated in said tabulating step. 



1 S 70. The method of Claim 1 1 , wherein said libraty of predetermined rotamer 

conformations is constructed ab initio, by a method comprising: computing portions of 
vibration-rotation potential energy surfaces of said side chain rotamers and deter mining, 
through exhaustive sampling, dihedral angles at which minima are found on said vibration- 
rotation potential energy surfaces. 

20 



71. The method of Claim 2, further comprising the step of storing all side chain 
rotamers having a calculated solution score that is lower than said threshold value in an 
array and eliminating a subset of side chain rotamers from said array using dead-end 
elimination where energy terms are partitioned according to; 



Template Energy + £ Intrinsic Energy* £ % Pairwise Energy 



residues residues residues 

J>9 



30 



72. The method of Claim 71» wherein dead-end elimination comprises elimination 
of a side chain rotamer r, of residue i, when the inequality 



35 

98 



BNSOOCID: <WO 01 16810A2_I_> 



WO 01/16810^ PCT/EP00/08504 



is true. 



73. The method of Claim 2, further comprising the step of eliminating a side chain 
5 rotamer pair using dead-end elimination, wherein a side chain rotamer pair comprises a first 

side chain rotamer corresponding to a first building block in said specified set of substitute 
building blocks and a second side chain rotamer corresponding to a second building block in 
said specified set of substitute building blocks. 

10 

74. The method of Claim 73, wherein dead-end elimination comprises eliminating a 
side chain rotamer pair that consists of rotamer r, of residue i, and rotamer 5. of residue j, 
when the inequalities: 



15 



and 

20 



25 



are true. 



75. The method of Claim 71, wherein the eliminating step is repeated until no 
additional side chain ro tamers are eliminated from said array by the process of dead-end 
30 elimination. 



76. The method of Claim 73, wherein the eliminating step is repeated until no 
additional side chain rotamer pairs are eliminated from said array by the process of dead-end 
35 elimination. 
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77. The method of Claim 2, wherein the calculated solution score additionally 

c mprises an entropy term that is determined by a difference between an entropy f a side 
chain rotamer in the solution when substituted into the atomic structure, and an entropy of 
the side chain rotamer in the solution when in a denatured state. 

5 

78. The method of Claim 77 wherein entropy contributions of said side chain 
rotamers are derived from experimental or empirical data. 

10 

79. The method of Claim 77, wherein the entropy term is computed using an 
iterative 

method. 

IS 

80. The method of Claim 79, wherein the entropy term is computed using mean 
field theory. 

20 

8 1 . The method of Claim 80, wherein the contribution of the rotamers to the entropy 
may be split into intrinsic and paiiwise terms. 

25 

82. The method of Claim 2 wherein said calculated solution score comprises a 
difference between a first value corresponding to said solution structure and a second value 
corresponding to a weighted average over a plurality of reference structures. 

30 

83. The method of Claim 1 wherein said target macromolecule consists of a 
plurality of structures at least one of which is represented by an atomic structure. 



35 84. A computer program product for use in conjunction with a computer, the 

computer program product comprising a computer readable storage medium and a computer 
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program mechanism embedded therein, the computer program mechanism comprising an 
optimizer module configured to choose a set of substitute building blocks for a set of 
positions in a target macromolecule according to whether a calculated solution score is 
lower than a threshold value, the computer program mechanism, upon receiving as input 
5 said set of positions: 

(a) specifying at least one substitute building block for each position in said set of 
positions to produce a specified set of substitute building blocks; 

(b) for each said substitute building block 

10 

i) determining at least one candidate conformer, 

ii) substituting coordinates of each said candidate conformer or portion 
thereof for coordinates of the building block or portion thereof at said position in an atomic 
structure of said target macromolecule; 

15 

(c) minimizing the value of a calculated energy term by adjusting the geometry of 
each said candidate conformer or portion thereof in order to obtain a solution structure; 

(d) calculating a solution score for said solution structure, wherein said solution 
2Q score comprises an entropic term; and 

(e) choosing said specified set of substitute building blocks if said calculated 
solution score is lower than a threshold value. 



85. The computer program product of Claim 84, wherein said macromolecule is a 
peptide or protein; said building blocks are amino acid residues; and each candidate 
conformer is a side chain rotamer selected from a plurality of side chain rotamers. 



30 

86. The computer program product of Claim 85 wherein said calculated solution 
score comprises a difference between a first value corresponding to said solution structure 
and a second value corresponding to a reference structure 

35 
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87. The computer program product of Claim 86, wherein said first value 
corresponding to said solution structure accounts for interactions between said side chain 

4 

retainer and said atomic structure, and a sum of interactions between all pairs of all possible 
side chain ro tamers. 

5 88. The computer program product of Claim 87, further comprising a step of 

rejecting a side chain rotamer when the value of said interactions between said side chain 
retainer and said atomic structure is greater than a threshold value, 

10 

89. The computer program product of Claim 86, wherein said reference structure is 
a denatured state of said solution structure. 

1 5 90. The computer program product of Claim 85, wherein the dihedral angles of said 

side chain rotamer are optimized in step (c). 

91 . The computer program product of Claim 85, wherein the positions of all main 
20 chain atoms of said atomic structure, and the positions of all atoms in amino acid side 
chains that are not included in said set of substitute building blocks are held fixed in said 
atomic structure. 



92. The computer program product of Claim 90, wherein the positions of all atoms 
in amino acid side chains on residues that are not at said set of positions are allowed to vary 
whilst the dihedral angles of said side chain rotamer are optimized. 

30 

93. The computer program product of Claim 90, wherein the positions of all main 
chain atoms of said atomic structure are allowed to vary whilst the dihedral angles of said 
side chain rotamer are optimized 

35 
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94. TTie computer program product of Claim 85, wherein said plurality of 
conformers is a library of predetermined side chain rotamer conformations. 



5 95. The computer program product of Claim 85, wherein at least one side chain 

conformer in said plurality of side chain retainers is derived from a continuous distribution 
of conformations. 



10 96. The computer program product of Claim 84, wherein said atomic structure 

includes a representation of the building blocks at each position in said set of positions; and 
said atomic structure was determined by a method selected from the group consisting of x- 
ray crystallography, nuclear magnetic resonance spectroscopy, electron microscopy, 
homology modeling, and ab initio molecular modeling. 



97. The computer program product of Claim 84, wherein said atomic structure is an 
X-ray crystal structure of a portion of said macromolecule that comprises said building 
blocks at each position. 



98. The computer program product of Claim 97, wherein said X-ray crystal 
structure was determined at a resolution of less than 4.0 Angstroms. 

25 

99. The computer program product of Claim 85, wherein said solution score is 
calculated using an empirical scoring function. 

30 

100. The computer program product of Claim 99, wherein said empirical scoring 
function is a sum of energy terms, comprising a template energy of said atomic structure 
held in a fixed geometry, an intrinsic energy of interaction between a side chain rotamer z 
said atomic structure held in a fixed geometry and a pairwise energy of interaction betwo 

35 possible pairs of side chain retainers in said substitute set of building blocks. 
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101. The computer program product of Claim 100, wherein the template energy of 
said atomic structure comprises at least one energy term selected from the group consisting 
of a molecular mechanics potential, a solvation energy, an empirical penalty function, and 
an entropic contribution. 

102, The computer program product of Claim 101, wherein the template energy of 
said atomic structure comprises a sum of terms whose coefficients are individually 
adjustable weighting factors. 



103, The computer program product of Claim 102, wherein said molecular 
mechanics potential comprises at least one energy term selected from the group consisting 
of bond length vibrations, bond angle bends, the hydrogen bond energy between pairs of 
1 5 hydrogen bond donor and acceptor atoms, an electrostatic interaction energy between pairs 
of charged atoms, and a van der Waals interaction energy between pairs of non-bonded 
atoms in said atomic structure. 



104, The computer program product of Claim 101, wherein said entropic 
contribution comprises at least one term selected from the group consisting of a main chain 
entropy term, a side chain rotation entropy term and a side chain vibration entropy term. 



25 



105. The computer program product of Claim 89, wherein said reference structure 
comprises said side chain rotamer embedded in an alanine based penta-peptide. 



30 106. The computer program product of Claim 89, wherein said reference structure 

comprises said side chain rotamer embedded in a fragment of protein taken from an atomic 
structure of a naturally occurring protein or an ensemble of fragments of protein, the 
populations of which are determined either from populations in the naturally occurring 
proteins or from computations establishing the potential energy of each fragment and 

35 integrating them into a partition function. 
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107. The computer program product of Claim 88, wherein the side chain retainers 
which are not rejected and that have a probability smaller than a predetermined probability 
threshold are rejected. 

5 

108. The computer program product of Claim 94, wherein said library of 
predetermined rotamer conformations is constructed by a method comprising: 

tabulating, for each of the twenty naturally occurring amino acids, a 
jq statistical distribution of observed amino acid side chain dihedral angles in a set of 

crystallographically determined protein structures; 

determining a Gaussian distribution of observed amino acid side chain 
dihedral angles tabulated in said tabulating step; and 

1 5 constructing amino acid side chain rotamers for each of the twenty naturally 

occurring amino acids using all combinations of Gaussian peaks around each amino 
acid side chain dihedral angle tabulated in said tabulating step. 

109. The computer program product of Claim 94, wherein said library of 
predetermined rotamer conformations is constructed ah initio, by a method comprising: 
computing portions of vibration-rotation potential energy surfaces of said side chain 
rotamers and determining, through exhaustive sampling* dihedral angles at which minima 
are found on said vibration-rotation potential energy surfaces. 

110. The computer program product of Claim 85, further comprising the step of 
storing all side chain rotamers having a solution score that is lower than said threshold value 
in an array and eliminating a subset of side chain rotamers from said array using dead-end 
elimination where energy terms are partitioned according to: 

Template Energy* £ Intrinsic Energy* £ £ Pairwise Energy 

residues residues residues 

. i i j>t 
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25 



30 
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111. The computer program product of Claim 1 1 0, wherein dead-end elimination 
comprises elimination of a side chain rotamer r, of residue i, when the inequality 



is true. 

1 12. The computer program product of Claim 85, further comprising the step of 
eliminating a side chain rotamer pair using dead-end elimination, wherein a side chain 
rotamer pair comprises a first side chain rotamer representing a first building block in said 
set of building blocks and a second side chain rotamer representing a second building block 
in said set of building blocks. 



113, The computer program product of Claim 112, wherein dead-end elimination 
20 comprises eliminating a side chain rotamer pair that consists of rotamer r. of residue i, and 
rotamer s, of residue j, when the inequalities: 

kmlj M.J 

25 



and 



30 



kmlj 



are true. 



35 
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1 14. The computer program product of Claim 1 10, wherein the eliminating step is 
repeated until no additi nal side chain rotamers are eliminated from said array by the 

pr cess of dead-end elimination. 

5 

115. The computer program product of Claim 1 12, wherein the eliminating step is 
repeated until no additional side chain rotamer pairs are eliminated from said array by the 
process of dead-end elimination. 

10 

1 16. The computer program product of Claim 85, wherein the solution score 
additionally comprises an entropy term that is determined by a difference between an 
entropy of a side chain rotamer in the solution when substituted into the atomic structure, 
and an entropy of the side chain rotamer in the solution when in a denatured state. 

15 

117. The computer program product of Claim 116, wherein entropy contributions of 
said side chain rotamers are derived from experimental or empirical data. 

20 

118. The computer program product of Claim 116, wherein the entropy term is 
computed using an iterative method. 

25 

119. The computer program product of Claim 118, wherein the entropy term is 
computed using mean field theory. 

30 

120. The computer program product of Claim 1 19, wherein the contribution of the 
rotamers to the entropy may be split into intrinsic and pairwise terms. 

35 121 . A system for choosing a set of substitute building blocks for a set f positions 

in a target macromoiecule according to whether a calculated solution score is 1 wer than a 
threshold value, comprising: 
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a central processing unit; 

an input device for inputting requests; 

an output device; 

5 

a memory; 

at least one bus connected to the central processing unit, the memory, the input 
device, and the output device; 

10 the memory storing an computer program mechanism comprising an optimizer 

module configured to choose the set of substitute building blocks, the computer program 
mechanism, upon receiving a request to choose the set of substitute building blocks, 

(a) specifying at least one substitute building block for each position in said set of 
j j positions to produce a specified set of substitute building blocks; 

(b) for each said substitute building block, 

i) determining at least one candidate confoimer; 

2 0 ii) substituting coordinates of each said candidate confonner or portion 

thereof for coordinates of the building block or portion thereof at said 
position in an atomic structure of said target macromolecule; 

(c) minimizing the value of a calculated energy term by adjusting the geometry of 
each said candidate confonner or portion thereof in order to obtain a solution structure; 

25 

(d) calculating a solution score for said solution structure, wherein said solution 
score comprises an entropic term; and 

(e) choosing said specified set of substitute building blocks if said calculated 
solution score is lower than a threshold value. 

30 

122. The system of Claim 121, wherein said macromolecule is a peptide or protein; 
said building blocks are amino acid residues; and each candidate confonner is a side chain 
35 ro tamer selected from a plurality of side chain rotamers. 
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123. The system of Claim 122 wherein said calculated solution score comprises a 
difference between a first value corresponding to said solution structure and a second value 
corresponding to a reference structure. 

5 

124. The system of Claim 123, wherein said first value corresponding to said 
solution structure accounts for interactions between said side chain ro tamer and said atomic 
structure, and a sum of interactions between all pairs of all possible side chain retainers. 

10 

125. The system of Claim 124, further comprising a step of rejecting a side chain 
rotamer when the value of said interactions between said side chain rotamer and said atomic 
structure is greater than a threshold value. 

15 

126. The system of Claim 123, wherein said reference structure is a denatured state 
of said solution structure. 



20 

127. The system of Claim 122, additionally comprising a step wherein the dihedral 
angles of said side chain rotamer are optimized in step (c). 



25 

^ 128. The system of Claim 122, wherein the positions of all main chain atoms of said 

atomic structure, and the positions of all atoms in amino acid side chains that are not 
included in said set of substitute building blocks are held fixed in said atomic 

structure. 

30 



129. Hie system of Claim 127, wherein the positions of all atoms in amino acid side 
chains on residues that are not at said set of positions are allowed to vary whilst the dihedral 
angles of said rotamer are optimized. 

35 
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130. The system of Claim 127, wherein the positions of all main chain atoms of said 
atomic structure are allowed to vary whilst the dihedral angles of said retainer are 
optimized 



131. The system of Claim 122, wherein said plurality of confonners is a library of 
predetermined rotamer conformations. 

10 132. The system of Claim 122, wherein at least one side chain rotamer in said 

plurality of side chain rotameis is derived from a continuous distribution of conformations. 



133. The system of Claim 121, wherein said atomic structure includes a 
representation of the building blocks at each position in said set of positions; and said 
atomic structure was determined by a method selected from the group consisting of x-ray 
crystallography, nuclear magnetic resonance spectroscopy, electron microscopy, homology 
modeling, and ab initio molecular modeling. 

20 

134. The system of Claim 121, wherein said atomic structure is an X-ray crystal 
structure of a portion of said macromolecule that comprises said building blocks at each 
position. 

25 

135. The system of Claim 124, wherein said X-ray crystal structure was determined 
at a resolution of less than 4.0 Angstroms. 

30 

136. The system of Claim 122, wherein said calculated solution score is obtained 
using an empirical scoring Amotion. 

35 

137. The system f Claim 136, wherein said empirical scoring function is a sum of 
energy terms, comprising a template energy of said atomic structure held in a fixed 
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geometry, an intrinsic energy of interaction between a side chain to tamer and said atomic 
structure held in a fixed geometry and a pairwise energy of interaction between possible 
pairs of side chain rotamers in said substitute set of building blocks. 

138. The system of Claim 137, wherein the template energy of said atomic structure 
comprises at least one energy term selected from the group consisting of a molecular 
mechanics potential, a solvation energy, an empirical penalty function, and an entropic 
contribution. 



139. The system of Claim 138, wherein the template energy of said atomic structure 
comprises a sum of terms whose coefficients are individually adjustable weighting factors. 

15 

140. The system of Claim 139, wherein said molecular mechanics potential 
comprises at least one energy term selected from the group consisting of bond length 
vibrations, bond angle bends, the hydrogen bond energy between pairs of hydrogen bond 

2Q donor and acceptor atoms, an electrostatic interaction energy between pairs of charged 
atoms, and a van der Waals interaction energy between pairs of non-bonded atoms in said 
atomic structure. 



25 141. The system of Claim 138, wherein said entropic contribution comprises at least 

one term selected from the group consisting of a main chain entropy term, a side chain 
rotation entropy term and a side chain vibration entropy term. 



142, The system of Claim 125, wherein said reference structure comprises said side 
chain rotamer substituted for a side chain rotamer in an alanine based penta-peptide. 



143. The system of Claim 125, wherein said reference structure comprises said side 
chain rotamer embedded in a fragment of protein taken from an atomic structure of a 
naturally occurring protein or an ensemble of fragments f protein, the p pulations of which 
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are determined either from populations in the naturally exurring proteins or from 
computations establishing the potential energy of each fragment and integrating them into a 
partition function. 

5 144. The system of Claim 124, wherein the side chain rotamers which are not 

rejected in and that have a probability smaller than a predetermined probability threshold are 
rejected. 

10 

145. The system of Claim 131. wherein said library of predetermined rotamer 
conformations is constructed by a method comprising: 

tabulating, for each of the twenty naturally occurring amino acids, a 
15 statistical distribution of observed amino acid side chain dihedral angles in a set of 

crystallographically determined protein structures; 

detennining a Gaussian distribution of observed amino acid side chain 
dihedral angles tabulated in said tabulating step; and 

20 constructing amino acid side chain rotamers for each of the twenty naturally 

occurring amino acids using all combinations of Gaussian peaks around each amino 
acid side chain dihedral angle tabulated in said tabulating step. 



25 146. The system of Claim 131, wherein said library of predetermined rotamer 

conformations is constructed ab initio, by a method comprising: computing portions of 
vibration-rotation potential energy surfaces of said side chain rotamers and deterniining, 
through exhaustive sampling, dihedral angles at which minima are found on said vibration- 
rotation potential energy surfaces. 



30 



147. The system of Claim 122, further comprising the step of storing all side chain 
rotamers having a solution score that is lower than said threshold value in an array and 
35 eliminating a subset of side chain rotamers from said array using dead-end eumination 
where energy terms are partitioned according to: 



112 



BNSDOCID: <WO 01 16610A2J_> 



WO 01 



PCT/EP00/08504 



Template Energy + £ Intrinsic Energy + £ £ Pairwise Energy 

residues residues rt&iduei 

i * j>i 

5 

148. The system of Claim 147, wherein dead-end elimination comprises elimination 
of a side chain retainer r, of residue i, when the inequality 



10 

-.template template y / paiiwe paiF>^c\ 

is true. 

15 

149. The system of Claim 122, further comprising the step of eliminating a side 
chain retainer pair using dead-end elimination, wherein a side chain ro tamer pair comprises 
a first side chain retainer representing a first building block in said specified set of substitute 
2Q building blocks and a second side chain retainer representing a second building block in 
said specified set of substitute building blocks. 



150. The system of Claim 149, wherein dead-end elimination comprises eliminating 
25 a side chain rotamcr pair that consists of retainer r, of residue i, and retainer s, of residue j, 
when the inequalities: 

+ Z > £ (w.) + I """.fawn* ) 
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are true. 



151. The system of Claim 147, wherein the eliminating step is repeated until no 
additional side chain retainers are eliminated from said array by the process of dead-end 
elimination. 



10 



15 



152. The system of Claim 149, wherein the eliminating step is repeated until no 
additional side chain rotamer pairs are eliminated from said array by the process of dead-end 
elimination. 



153. The system of Claim 122. wherein the solution score additionally comprises an 
entropy term that is determined by a difference between an entropy of a side chain rotamer 
in the solution when substituted into the atomic structure, and an entropy of the side chain 
rotamer in the solution when in a denatured state. 



20 



154. The system of Claim 153, wherein entropy contributions of said side chain 
rotamers are derived from experimental or empirical data. 

25 

1 55. The system of Claim 153, wherein Ihe entropy term is computed using an 
iterative method 

30 

156. The system of Claim 153, wherein the entropy term is computed using mean 
field theory. 



157. The system of Claim 154, wherein the contribution of the rotamere to the 
entropy may be split into intrinsic and pairwise terms. 
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158. The method of Claim 1 additionally comprising the step of synthesizing at 
least one of said solution structures which has a calculated solution score lower than said 
threshold value. 

5 

159. The method of Claim 158 additionally comprising the step of screening each of 
said solution structures that has been synthesized against an assay to test for activity. 

*® 160. The method of Claim 1 wherein at least one of said building blocks in said set 

of building blocks contacts a binding partner of said target raacromoleculc when bound to 
said target macromolecule. 



161. The method of Claim 2 wherein said calculated solution score comprises a 
difference between a first value corresponding to said solution structure and a second value 
corresponding to a denatured state of said solution structure, 

wherein said calculated solution score is calculated using an empirical scoring 
2 ^ function that comprises a sum of energy terms, including an energy of said atomic structure 
held in a fixed geometry, an intrinsic energy of interaction between a candidate side chain 
rotamer and said atomic structure held in a fixed geometry and a pairwise energy of 
interaction between possible pairs of side chain rotamers in said substitute set of building 
blocks, 

25 

wherein any one of said energy of said atomic structure, said intrinsic energy of 
interaction and said pairwise energy of interaction comprises at least one energy term 
selected from the group consisting of a molecular mechanics potential, a solvation energy, 
an empirical penalty function, and an entropic contribution, and 

30 

wherein said entropic contribution is computed using mean field theory. 



35 
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SEQUENCE LISTING 

<110> The European Molecular Biology Laboratory (EMBL) 

<120> A COMPUTER-BASED METHOD FOR 

MACROMOLECULAR ENGINEERING AND DESIGN 



<130> 9882-017-22B 

<150> US 09/387,741 
<1S1> 1999-08-31 

<160> 20 

<170> FaetSEQ for Windows Version 4.0 

<210> 1 
<211> 6 
<212> PRT 

*213> Artificial Sequence 

<22Q> J _ . 

<223> Conceptual sequence for computer modeling 

c400> 1 

Ala Val He Leu Phe Trp 
1 5 



<210> 2 
<21X> 14 
<212> PRT 

<213> Artificial Sequence 
c220> 

<223> Conceptual sequence for computer modeling 
<400> 2 

Ala Val He Leu Gly Asp Aen Ser Thr Glu Lys Arg Tyr Trp 



<210> 3 
<211> 4 
<212> PRT 

<213> Artificial Sequence 
<220> 

<223> Conceptual sequence 
<400> 3 

Ala Val He Leu 



for computer modeling 



<210> 4 
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<211> 5 
«212> PRT 

<213> Artificial Sequence 
<220> 

<223> Conceptual sequence for computer modeling 
<400> 4 

Ala Val He Leu Phe 

1 5 



c210> 5 
c211> 4 
<212> PRT 

<213> Artificial Sequence 
<220> 

<223> Conceptual sequence for computer modeling 
<400> s 

Ala Val Leu Phe 
1 



<210> 6 
<211» 8 
<212> PRT 

<213> Artificial Sequence 
<220> 

<223> Conceptual sequence for computer modeling 

<400> 6 

Val Ala val Met Leu Leu Trp Val 
1 5 



<210> 7 
<211> 9 
<212> PRT 

<213> Artificial Sequence 
<220> 

<223> Conceptual sequence for computer modeling 
<400> 7 

He Val lie He Leu Leu Val He Val 
1 5 



<210> 8 
<211> 4 
<212> PRT 

<213s» Artificial Sequence 

<220> 

<223> Conceptual sequence for computer modeling 

2 
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<400> 8 

Val Gly Ser Lye 

1 



<2XO> 9 
<2X1> 9 
<212> PRT 

<2X3> Artificial Sequence 
<220> 

<223> Conceptual sequence for computer modeling 
c400> 9 

Leu Val lie Val Leu Leu Val lie Val 
1 5 



<210> 10 
<211> 9 
<212> PRT 

<2X3? Artificial Sequence 
*220> 

<223> Conceptual sequence for computer modeling 
<400> 10 

Val Val He lie Leu Leu Val He Val 
1 S 



<210> 11 
<211> 9 
<212> PRT 

<2X3;> Artificial Sequence 
<220> 

<223> Conceptual sequence for computer modeling • 
<400> 11 

He Val He He Leu Leu Val val Val 
1 5 



<210> 12 
<211> 9 
<212> PRT 

<213> Artificial Sequence 
<220> 

<223> Conceptual sequence for computer modeling 

<400> 12 

Leu He He Val Leu Leu Val He val 
X 5 



<210> X3 



3 
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<211> 9 
<212> PRT 

<213> Artificial Sequence 
<220> 

<223> Conceptual eequence for computer modeling 
<400> 13 

He Val Val He Leu Leu Val He Val 
1 S 



<2\0> 14 
<211> 9 
<212> PRT 

<213> Artificial Sequence 
<220> 

<223> Conceptual eequence for computer modeling 
<400> 14 
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